C++ TTree::GetEntries方法代码示例

void macroEmuTree() {
  string inputline;
  string outputline;
  string blankline;
  ifstream myfile ("listofsamples.txt");
  if (myfile.is_open())
    {
      while ( myfile.good() )
	{
	  getline (myfile,inputline);
	  getline (myfile,outputline);
	  getline (myfile,blankline);

	  //gSystem->Load("$ROOTSYS/test/libEvent");
	  cout << "File : " << inputline << endl;
	  //Get old file, old tree and set top branch address
	  TString inputfilepath =inputline; 
	  TString outputfilepath = outputline; 
	  //TDCacheFile *oldfile = new TDCacheFile(inputfilepath);
	  TFile *oldfile = new TFile(inputfilepath);
	  TTree *oldtree = (TTree*)oldfile->Get("gsfcheckerjob/tree");
	  Long64_t nentries = oldtree->GetEntries();

          // Declaration of leaf types
          UInt_t          runnumber;
          UInt_t          eventnumber;
          UInt_t          luminosityBlock;
          Int_t           HLT_Mu22_Photon22_CaloIdL;
          Int_t           HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;
          Int_t           HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;
          Int_t           prescale_HLT_Mu22_Photon22_CaloIdL;
          Int_t           prescale_HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;
          Int_t           prescale_HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;
          Float_t         rho;
          Float_t         pfmet;
          Int_t           pvsize;
          Int_t           JetColl_size;
          Float_t         Jet_pt[100];   //[JetColl_size]
          Int_t           muon_size;
          Float_t         muon_pt[100];   //[muon_size]
          Float_t         muon_ptError[100];   //[muon_size]
          Float_t         muon_eta[100];   //[muon_size]
          Float_t         muon_phi[100];   //[muon_size]
          Int_t           muon_charge[100];   //[muon_size]
          Int_t           muon_nhitspixel[100];   //[muon_size]
          Int_t           muon_nhitstrack[100];   //[muon_size]
          Int_t           muon_nhitsmuons[100];   //[muon_size]
          Int_t           muon_nlayerswithhits[100];   //[muon_size]
          Int_t           muon_nSegmentMatch[100];   //[muon_size]
          Bool_t          muon_isTrackerMuon[100];   //[muon_size]
          Float_t         muon_normChi2[100];   //[muon_size]
          Float_t         muon_dz_beamSpot[100];   //[muon_size]
          Float_t         muon_dz_firstPVtx[100];   //[muon_size]
          Float_t         muon_dxy_cmsCenter[100];   //[muon_size]
          Float_t         muon_dxy_beamSpot[100];   //[muon_size]
          Float_t         muon_dxy_firstPVtx[100];   //[muon_size]
          Float_t         muon_trackIso03[100];   //[muon_size]
          Float_t         muon_emIso03[100];   //[muon_size]
          Float_t         muon_hadIso03[100];   //[muon_size]
          Int_t           gsf_size;
          Float_t         gsf_eta[100];   //[gsf_size]
          Float_t         gsf_phi[100];   //[gsf_size]
          Float_t         gsf_theta[100];   //[gsf_size]
          Int_t           gsf_charge[100];   //[gsf_size]
          Float_t         gsf_sigmaetaeta[100];   //[gsf_size]
          Float_t         gsf_sigmaIetaIeta[100];   //[gsf_size]
          Float_t         gsf_dxy_firstPVtx[100];   //[gsf_size]
          Float_t         gsf_dz_beamSpot[100];   //[gsf_size]
          Float_t         gsf_dz_firstPVtx[100];   //[gsf_size]
          Int_t           gsf_nLostInnerHits[100];   //[gsf_size]
          Float_t         gsf_deltaeta[100];   //[gsf_size]
          Float_t         gsf_deltaphi[100];   //[gsf_size]
          Float_t         gsf_hovere[100];   //[gsf_size]
          Float_t         gsf_trackiso[100];   //[gsf_size]
          Float_t         gsf_ecaliso[100];   //[gsf_size]
          Float_t         gsf_hcaliso1[100];   //[gsf_size]
          Float_t         gsf_hcaliso2[100];   //[gsf_size]
          Bool_t          gsf_isecaldriven[100];   //[gsf_size]
          Float_t         gsfsc_e[100];   //[gsf_size]
          Float_t         gsfsc_eta[100];   //[gsf_size]
          Float_t         gsfsc_phi[100];   //[gsf_size]
          Float_t         gsf_e2x5overe5x5[100];   //[gsf_size]
          Float_t         gsf_e1x5overe5x5[100];   //[gsf_size]
          Float_t         gsf_gsfet[100];   //[gsf_size]
          Float_t         genPair_mass;
          Int_t           trueNVtx;
       
          // List of branches
          TBranch        *b_runnumber;   //!
          TBranch        *b_eventnumber;   //!
          TBranch        *b_luminosityBlock;   //!
          TBranch        *b_HLT_Mu22_Photon22_CaloIdL;   //!
          TBranch        *b_HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;   //!
          TBranch        *b_HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;   //!
          TBranch        *b_prescale_HLT_Mu22_Photon22_CaloIdL;   //!
          TBranch        *b_prescale_HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;   //!
          TBranch        *b_prescale_HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL;   //!
          TBranch        *b_rho;   //!
          TBranch        *b_pfmet;   //!
          TBranch        *b_pvsize;   //!
//.........这里部分代码省略.........

// Simple example of reading a generated Root file
void testgeo(char *filename=NULL)
{

  // Clear global scope
  gROOT->Reset();

  // Load the library with class dictionary info
  // (create with "gmake shared")
  char* wcsimdirenv;
  wcsimdirenv = getenv ("WCSIMDIR");
  if(wcsimdirenv !=  NULL){
    gSystem->Load("${WCSIMDIR}/libWCSimRoot.so");
  }else{
    gSystem->Load("../libWCSimRoot.so");
  }

  // Open the file
  TFile *file; 
  if (filename==NULL){
    file = new TFile("../wcsim.root","read");
  }else{
    file = new TFile(filename,"read");
  }
  if (!file->IsOpen()){
    cout << "Error, could not open input file: " << filename << endl;
    return -1;
  }
  
  // Get the a pointer to the tree from the file
  TTree *gtree = (TTree*)file->Get("wcsimGeoT");
  
  // Get the number of events
  int nevent = gtree->GetEntries();
  printf("geo nevent %d\n",nevent);
  
  // Create a WCSimRootGeom to put stuff from the tree in

  WCSimRootGeom* wcsimrootgeom = new WCSimRootGeom();

  // Set the branch address for reading from the tree
  TBranch *branch = gtree->GetBranch("wcsimrootgeom");
  branch->SetAddress(&wcsimrootgeom);

  // Now loop over "events"  (should be only one for geo tree)
  int ev;
  for (ev=0;ev<nevent; ev++)  {
      // Read the event from the tree into the WCSimRootGeom instance
      gtree->GetEntry(ev);
      printf("Cyl radius %f\n", wcsimrootgeom->GetWCCylRadius());
      printf("Cyl length %f\n", wcsimrootgeom->GetWCCylLength());
      printf("PMT radius %f\n", wcsimrootgeom->GetWCPMTRadius());
      printf("Offset x y z %f %f %f\n", wcsimrootgeom->GetWCOffset(0),
	     wcsimrootgeom->GetWCOffset(1),wcsimrootgeom->GetWCOffset(2));
      int numpmt = wcsimrootgeom->GetWCNumPMT();
      printf("Num PMTs %d\n", numpmt);

      int i;
      for (i=0;i<((numpmt<20)?numpmt:20);i++){
	WCSimRootPMT pmt;
	pmt = wcsimrootgeom->GetPMT(i);
	printf ("pmt %d %d %d\n",i,pmt.GetTubeNo(), pmt.GetCylLoc());
	printf ("position: %f %f %f\n", pmt.GetPosition(0),
		pmt.GetPosition(1),pmt.GetPosition(2));
	printf ("orientation: %f %f %f\n", pmt.GetOrientation(0),
		pmt.GetOrientation(1),pmt.GetOrientation(2));
      }
            
    } // End of loop over events

}

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

  }
    
  //-- initialize all variables
  for(int icent = 0; icent<NCENT; icent++){
    for(int iEP = 0; iEP < NEPSymm; iEP++){
      if( iEP != EPSymmBin ) continue;
      for(int iqn = 0; iqn < NQN; iqn++){

	VNDetX_0[icent][iEP][iqn]     = 0.;
	VNDetY_0[icent][iEP][iqn]     = 0.;
	VNDetX_1[icent][iEP][iqn]     = 0.;
	VNDetY_1[icent][iEP][iqn]     = 0.;
	VNDetX_full[icent][iEP][iqn]  = 0.;
	VNDetY_full[icent][iEP][iqn]  = 0.;

	evtMult_0[icent][iEP][iqn]    = 0;
	evtMult_1[icent][iEP][iqn]    = 0;
	evtMult_full[icent][iEP][iqn] = 0;

	Nevents[icent][iEP][iqn]      = 0;
	NFails[icent][iEP][iqn]       = 0;

      }
    }
  }
    
    
  //
  // Calculate Vn_det
  //
    
  cout<<"Begin DETECTOR loop, contains "<<tree->GetEntries()<<" Events"<<endl;
    
  int N;
  if(testrun) N = 10000; 
  else        N = tree->GetEntries();


  //-- Begin event loop
  for(int ievent = 0; ievent < N; ievent++) {
        
    if((ievent+1)% 500000 == 0) cout << "Processing Event " << ievent+1 << "\t" << (100.*(ievent+1)/N) << "% Completed" << endl;
        
    tree->GetEntry(ievent);
      
    //-- Vertex Cut
    if(TMath::Abs(vtx) < 3. || TMath::Abs(vtx) > 15.) continue;
      
    //-- Calculate centbin
    if( centval > cent_max[NCENT-1]) continue;
    int icent = hCentBins.FindBin(centval)-1;

    //-- begin EP loop
    for(int iEP = 0; iEP < NEP; iEP++){

      int EPbin  = EPSymmPartnerBin[iEP];
      if( EPbin != EPSymmBin ) continue;

      //-- Calculate qbin
      double qx    = qnHFx_EP[iEP];
      double qy    = qnHFy_EP[iEP];
      double qxDet = hqnHFDet_x[iEP]->GetBinContent(icent+1);
      double qyDet = hqnHFDet_y[iEP]->GetBinContent(icent+1);
      double sumET = sumET_EP[iEP];

void run(std::string inputFile, int max_iterations, int entries_per_kernel = 100) {
  
  TFile *file2 = new TFile(inputFile.c_str());
  
 
  std::vector<double>* samples = new std::vector<double>;
  double amplitudeTruth;
  TTree *tree = (TTree*)file2->Get("Samples");
  tree->SetBranchAddress("amplitudeTruth",      &amplitudeTruth);
  tree->SetBranchAddress("samples",             &samples);
  int nentries = tree->GetEntries();
  
  float time_shift = 13. ; //---- default is 13
  float pedestal_shift = 0.;
  
  
  float return_chi2 = -99;
  float best_pedestal = 0;
  float best_chi2 = 0;
  
  
  std::vector<TH1F*> v_pulses;
  std::vector<TH1F*> v_amplitudes_reco;
  
  fout->cd();
  TTree* newtree = (TTree*) tree->CloneTree(0); //("RecoAndSim");
  newtree->SetName("RecoAndSim");
  
  std::vector <double> samplesReco;
  std::vector < std::vector<double> > complete_samplesReco;
  std::vector <double> complete_chi2;
  std::vector <double> complete_pedestal;
  
  int ipulseintime = 0;
  newtree->Branch("chi2",   &return_chi2, "chi2/F");
  newtree->Branch("samplesReco",   &samplesReco);
  newtree->Branch("ipulseintime",  ipulseintime,  "ipulseintime/I");
  
  newtree->Branch("complete_samplesReco",   &complete_samplesReco);
  newtree->Branch("complete_chi2",          &complete_chi2);
  newtree->Branch("complete_pedestal",          &complete_pedestal);
  newtree->Branch("best_pedestal",   &best_pedestal, "best_pedestal/F");
  newtree->Branch("best_chi2",   &best_chi2, "best_chi2/F");
  
  
  int totalNumberOfBxActive = 10;
  
  for (unsigned int ibx=0; ibx<totalNumberOfBxActive; ++ibx) {
    samplesReco.push_back(0.);
  }
  
  
  v_amplitudes_reco.clear();
  
  
  struct Args {
      SampleVector samples;
      SampleMatrix samplecor;
      double pederr;
      BXVector bxs;
      FullSampleVector fullpulse;
      FullSampleMatrix fullpulsecov;

      Args(SampleVector const& samples,
           SampleMatrix const& samplecor,
           double pederr,
           BXVector const& bxs,
           FullSampleVector fullpulse,
           FullSampleMatrix fullpulsecov)
        : samples(samples), samplecor(samplecor), pederr(pederr), bxs(bxs), 
          fullpulse(fullpulse), fullpulsecov(fullpulsecov)
      {}
  };

  struct Output {
      double chi2;
      double ampl;
      int status;

      Output(double chi2, double ampl, int status) 
          : chi2{chi2}, ampl{ampl}, status{status}
      {}
  };
  
  std::cout 
            << "max_iterations: " << max_iterations << std::endl
            << "entries_per_kernel: " << entries_per_kernel << std::endl;
  
  for (auto it=0; it<max_iterations; ++it) {
      // vector of input parameters to the kernel
      std::vector<Args> vargs;

      for (int ie=0; ie<entries_per_kernel; ++ie) {
          tree->GetEntry(ie % tree->GetEntries());
          for (int i=0; i<NSAMPLES; ++i)
              amplitudes[i] = samples->at(i);

          double pedval = 0.;
          double pedrms = 1.0;
          vargs.emplace_back(amplitudes, noisecor, pedrms, activeBX, fullpulse, fullpulsecov);
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  ApplyLHCbStyle();
  double maxrate;
  switch(argc)
  {
    case 2:
      maxrate=-1;
      break;
    case 3:
      if(atoi(argv[2]))
        maxrate=atoi(argv[2]);
      else
      {
        cout << "Argument not a number" << endl;
        return 1;
      }
      break;
    default:
      cout << "Usage: " << argv[0] << " <Filename> [<Max rate>]" << endl;
      return 1;
  }
  if(system("[ -a figs ]") != 0)
  {
    cout << "Directory ./figs/ does not exist" << endl;
    const int mkdirerr = system("mkdir figs");
    if(mkdirerr != 0)
    {
      cout << "Error code " << mkdirerr << " while creating ./figs/" << endl;
    }
    else
    {
      cout << "Directory ./figs/ created" << endl;
    }
  }
  TFile* file = new TFile(argv[1]);
  if(!file->IsOpen())
  {
    cout << "Exiting" << endl;
    return 1;
  }
  TTree* tree = (TTree*)file->Get("metatree");
  int iHV; tree->SetBranchAddress("HV", &iHV);
  vector<int> HVs;
  string HV;
  int n = tree->GetEntries();
  TCanvas* mapscan, * histcan, * graphcan, * indivgraphcan[4];
  TPad* mapspad[4], * histpad[4];
  TH1D* hist;
  TH2D* map;
  string title[4] = {"JB", "JT", "ST", "SB"};
//  string title[4] = {"FA0026", "DA0024", "FA0019", "FA0006"};
  TGraphErrors* graph[4];
  TMultiGraph* mg;
  TLegend* lg;
  float xlo, xhi, ylo, yhi;
//  rainbowgradient();
  heatmapgradient();
//  DarkBodyRadiator();
//  DarkBodyRadiator2();
//  redbluegradient();
//  heatmapRB();
  string plotname;
  stringstream plotnamestream;
  double* x[4], * y[4], * xe[4], * ye[4];
  for(int j = 0; j < 4; j++)
  {
    x[j]  = new double[n];
    y[j]  = new double[n];
    xe[j] = new double[n];
    ye[j] = new double[n];
  }

  Double_t maxima[4];
  cout << "Entering loop of " << n << " events." << endl;
  for(int i = 0; i < n; i++)
  {
    tree->GetEntry(i);
    HVs.push_back(iHV);
    HV = itoa(iHV);
    for(int j = 0; j < 4; j++)
    {
      map = (TH2D*)file->Get(("DarkMap"+HV+"V"+itoa(j)).c_str());
      maxima[j] = map->GetMaximum();
    }
  }
  sort(HVs.begin(), HVs.end());
  if(maxrate==-1)
  {
    for(int i = 1; i < 4; i++)
    {
      maxrate = TMath::Max(maxima[i-1], maxima[i]);
      maxima[i]=maxrate;
    }
  }
//  cout << "Plotting with a maximum rate of " << maxrate << endl;
  for(int i = 0; i < n; i++)
  {
    HV = itoa(HVs[i]);
    plotnamestream << "DarkRate" << setfill('0') << setw(4) << HVs[i];
//.........这里部分代码省略.........

//void Plot_Yield(){
int main(){
	//Here are where you need to change
	//const TString Com = "DpTh40";
	TString Target = "He3"; cerr<<"--- Target (H2,He3,He4,C12,Ca40 and Ca48) = "; cin >> Target;
	TString Kin = "Kin3.1"; cerr<<"--- Kin (e.g. 3.1) = "; cin >> Kin; Kin="Kin"+Kin;
	TString Arm = "L";      cerr<<"--- Arm (L or R) = "; cin >> Arm;
	TString Com = "J02";    cerr<<"--- Comment = "; cin >> Com;

	/*Read In RunNo{{{*/
	TString sTarget="", sKin = "";
	if(Target.Contains("H")){
	    sTarget = Target; sTarget.ReplaceAll("H","h");
	}
   	if(Target.Contains("C")){
	    sTarget = Target; sTarget.ReplaceAll("C","c");
	} 
	sKin=Kin; sKin.ReplaceAll("Kin","kin");
	
	TString Input_Parameter = Form("%s/xs_input/%sHRS/xs_%s_%s.dat",XGT2_DIR.Data(),Arm.Data(),sTarget.Data(),sKin.Data());
	ifstream input_par(Input_Parameter.Data());
	TString one_line;
	
	int RunNo_Temp = gGet_InputFile_Var(Input_Parameter,0).Atoi();
    vector<int> vRunNoChain = gGet_RunNoChain(RunNo_Temp, gGet_InputFile_Var(Input_Parameter,15));  

	const int N=vRunNoChain.size();
    int RunNo[N];
	for(int i=0;i<N;i++)
		RunNo[i] = vRunNoChain[i];
	/*}}}*/

	TString InputFileName = Form("%s/results/%s/%s/%s_%s_%s_Yield_%s.root",
			XGT2_DIR.Data(),Target.Data(),Com.Data(),Arm.Data(),Kin.Data(),Target.Data(),Com.Data());
	cerr<<" &&& Reading in RootFile --> "<<InputFileName.Data()<<endl;

	double Ep,xbj;
	double Y_EX[N], Y_EX_Stat[N], Y_EX_Sys[N], Y_MC[N],Y_MC_Stat[N],Y_MC_Sys[N]; 

	TFile *f1 = new TFile(InputFileName.Data(),"r");
	if(f1->IsZombie()){
	   cerr<<"*** I cannot fine ROOTFile: "<<InputFileName.Data()<<" ***"<<endl;
	   exit(-1);
	}
		
	TTree *T = (TTree*) gDirectory->Get("T");

	const int N_Evt = T->GetEntries();

	T->SetBranchAddress("Ep", &Ep);
	T->SetBranchAddress("xbj", &xbj);
	for(int i=0;i<N;i++){
		T->SetBranchAddress(Form("Y_EX_%d",     RunNo[i]), &Y_EX[i]);
		T->SetBranchAddress(Form("Y_EX_Stat_%d",RunNo[i]), &Y_EX_Stat[i]);
		T->SetBranchAddress(Form("Y_EX_Sys_%d" ,RunNo[i]), &Y_EX_Sys[i]);
		T->SetBranchAddress(Form("Y_MC_%d",     RunNo[i]), &Y_MC[i]);
		T->SetBranchAddress(Form("Y_MC_Stat_%d",RunNo[i]), &Y_MC_Stat[i]);
		T->SetBranchAddress(Form("Y_MC_Sys_%d" ,RunNo[i]), &Y_MC_Sys[i]);
	}

	double Y_R[N][N_Evt], Y_R_Err[N][N_Evt];
	double Ep_Bin[N_Evt],xbj_Bin[N_Evt];
	double Ep_Max = -1000.0, Ep_Min = 1000.0, xbj_Max = -1000.0, xbj_Min = 1000.0, Y_EX_Err=0.0, Y_MC_Err = 0.0;
	for(int j=0;j<N_Evt;j++){
		T->GetEntry(j);
		Ep_Bin[j] = Ep;
		if(Ep>Ep_Max) Ep_Max = Ep;
		if(Ep<Ep_Min) Ep_Min = Ep;

		xbj_Bin[j] = xbj;
		if(xbj>xbj_Max) xbj_Max = xbj;
		if(xbj<xbj_Min) xbj_Min = xbj;


		for(int i=0;i<N;i++){
			//cerr<<Form("Run#%d: Ep = %f, Y_R = %e (%e)", RunNo[i],Ep_Bin[j], Y_EX[i],Y_EX_Stat[i])<<endl;
			Y_EX_Err = sqrt( pow(Y_EX_Stat[i],2)+pow(Y_EX_Sys[i],2) ); 
			Y_MC_Err = sqrt( pow(Y_MC_Stat[i],2)+pow(Y_MC_Sys[i],2) );
			if(Y_EX[i]>1e-30&&Y_MC[i]>1e-30){
				Y_R[i][j] = Y_EX[i] / Y_MC[i];
				Y_R_Err[i][j] = Y_R[i][j]*sqrt( pow(Y_EX_Err/Y_EX[i],2)+pow(Y_MC_Err/Y_MC[i],2) );
			}
			else{
				Y_R[i][j] = 0.0;
				Y_R_Err[i][j] = 0.0;
			}
		}
	}

	/*Yield  Ratio{{{*/
	double Y_Min =0, Y_Max = 2.0;
    if(N>10)
		Y_Min = -1.0;
	TCanvas *c3 = new TCanvas("c3","c3",1200,800);
	TH2F* r = new TH2F("r","",1000,xbj_Min-0.2,xbj_Max+0.2,1000,Y_Min, Y_Max);
	r->SetStats(kFALSE);
	r->SetXTitle("X");
	r->GetXaxis()->CenterTitle(1);
	r->GetXaxis()->SetTitleFont(32);
	r->GetXaxis()->SetTitleSize(0.06);
//.........这里部分代码省略.........

int main(int argc,char** argv){

    int pserr = getargs(argc,argv);

    TString Target = gTarget;
    TString aKin = gKin; TString Kin = Form("Kin%s",aKin.Data());
    TString Arm = gArm;
    TString Com = gCom;


    /*    //Here are where you need to change*/
    ////const TString Com = "DpTh40";
    //TString Target = "He3"; cerr<<"--- Target (H2,He3,He4,C12,Ca40 and Ca48) = "; cin >> Target;
    //TString Kin = "Kin3.1"; cerr<<"--- Kin (e.g. 3.1) = "; cin >> Kin; Kin="Kin"+Kin;
    //TString Arm = "L";      cerr<<"--- Arm (L or R) = "; cin >> Arm;
    //TString Com = "J02";    cerr<<"--- Comment = "; cin >> Com;

    if(Arm!="L"&&Arm!="R") exit(-3);

    /*Read In RunNo{{{*/
    TString sTarget="", sKin = "";
    if(Target.Contains("H")){
        sTarget = Target; sTarget.ReplaceAll("H","h");
    }
    if(Target.Contains("C")){
        sTarget = Target; sTarget.ReplaceAll("C","c");
    } 
    sKin=Kin; sKin.ReplaceAll("Kin","kin");

    //TString Input_Parameter = Form("%s/xs_input/%sHRS/xs_%s_%s.dat",XGT2_DIR.Data(),Arm.Data(),sTarget.Data(),sKin.Data());
    TString Input_Parameter = Form("%s/xs_input/%sHRS_FULL/xs_%s_%s.dat",XGT2_DIR.Data(),Arm.Data(),sTarget.Data(),sKin.Data());
    ifstream input_par(Input_Parameter.Data());
    TString one_line;

    int RunNo_Temp = gGet_InputFile_Var(Input_Parameter,0).Atoi();
    vector<int> vRunNoChain = gGet_RunNoChain(RunNo_Temp, gGet_InputFile_Var(Input_Parameter,15));  

    const int N_Run=vRunNoChain.size();
    int RunNo[N_Run];
    cout<<endl<<"--- RunList =";
    for(int i=0;i<N_Run;i++){
        RunNo[i] = vRunNoChain[i];
        cout<<"  "<<RunNo[i];
    }
    cout<<endl;
    /*}}}*/

    TString InputFilePath = Form("%sresults/%s/%s/", XGT2_DIR.Data(),Target.Data(),Com.Data());

    TString InputFileName = Form("%s/results/%s/%s/%s_%s_%s_Yield_%s.root",
            XGT2_DIR.Data(),Target.Data(),Com.Data(),Arm.Data(),Kin.Data(),Target.Data(),Com.Data());
    if(!(gSystem->FindFile(InputFilePath.Data(), InputFileName))){
        cerr<<" *** RootFile does exist--> "<<InputFileName.Data()<<endl;
        return -1;
    }

    cerr<<" &&& Reading in RootFile --> "<<InputFileName.Data()<<endl;

    double Ep,xbj;
    double XS_EX[N_Run], XS_MC[N_Run], XS_EX_Stat[N_Run], XS_EX_Sys[N_Run], Xbj_EX[N_Run];

   TString InputDataName = Form("%s/results/%s/%s/%s_%s_%s_Yield.dat_%s",/*{{{*/
            XGT2_DIR.Data(),Target.Data(),Com.Data(),Arm.Data(),Kin.Data(),Target.Data(),Com.Data());
  
    if(!(gSystem->FindFile(InputFilePath.Data(), InputDataName))){
        cerr<<" *** RootFile does exist--> "<<InputDataName.Data()<<endl;
        return -2;
    }
   cerr<<" &&& Reading in DataFile --> "<<InputDataName.Data()<<endl;
   ifstream run_1(InputDataName.Data());

    double dummy;
    TString com;
    com.ReadLine(run_1);
    vector<double> vXS_MC, vXbj;
    double xs_mc, xbj_mc;
    while ( run_1>>dummy>>dummy>>dummy>>xs_mc >> xbj_mc
                 >>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy
                 >>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy>>dummy){

        cout<<Form("--- From Data: xbj=%f,  XS_MC=%e", xbj_mc, xs_mc)<<endl;

        vXS_MC.push_back(xs_mc);
        vXbj.push_back(xbj_mc);
    }/*}}}*/

    TFile *f1 = new TFile(InputFileName.Data(),"r");
    if(f1->IsZombie()){
        cerr<<"*** I cannot fine ROOTFile: "<<InputFileName.Data()<<" ***"<<endl;
        exit(-1);
    }
    TTree *T = (TTree*) gDirectory->Get("T");
    const int N_Evt = T->GetEntries();
    T->SetBranchAddress("Ep", &Ep);
    T->SetBranchAddress("xbj", &xbj);
    for(int i=0;i<N_Run;i++){
        cout<<Form("--- For Run#%d", vRunNoChain[i])<<endl;
        T->SetBranchAddress(Form("XS_%d",     vRunNoChain[i]), &XS_EX[i]);
        T->SetBranchAddress(Form("XS_Stat_%d",vRunNoChain[i]), &XS_EX_Stat[i]);
        T->SetBranchAddress(Form("XS_Sys_%d" ,vRunNoChain[i]), &XS_EX_Sys[i]);
//.........这里部分代码省略.........

void listMismatchedEvents(TString inputFile = "L1UnpackedPureEmulator.root", testObject type = Jets)
{
  TFile *inFile = TFile::Open(inputFile);
  TTree *emulatorResults = (TTree*)inFile->Get("EmulatorResults/L1UpgradeTree");
  TTree *unpackerResults = (TTree*)inFile->Get("UnpackerResults/L1UpgradeTree");

  int event, run, lumi;

  emulatorResults->SetBranchAddress("event",&event);
  emulatorResults->SetBranchAddress("run",&run);
  emulatorResults->SetBranchAddress("lumi",&lumi);

  std::vector<int> *e_hwPt =0;
  std::vector<int> *e_hwEta =0;
  std::vector<int> *e_hwPhi =0;
  std::vector<int> *e_bx =0;
  int e_N;
  std::vector<int> *u_hwPt =0;
  std::vector<int> *u_hwEta =0;
  std::vector<int> *u_hwPhi =0;
  std::vector<int> *u_bx =0;
  int u_N;

  TString prefix;
  TString nPrefix;
  switch(type) {
  case Jets:
    prefix = "jet";
    nPrefix = "Jet";
    break;
  case EGammas:
    prefix = "egamma";
    nPrefix = "Egamma";
    break;
  case Taus:
    prefix = "tau";
    nPrefix = "Tau";
    break;
  case Centrality:
    prefix = "hfring";
    nPrefix = "Hfring";
    break;
  default:
    prefix = "jet";
    nPrefix = "Jet";
    break;
  }

  emulatorResults->SetBranchAddress(prefix + "_hwPt",&e_hwPt);
  emulatorResults->SetBranchAddress(prefix + "_hwEta",&e_hwEta);
  emulatorResults->SetBranchAddress(prefix + "_hwPhi",&e_hwPhi);
  emulatorResults->SetBranchAddress(prefix + "_bx",&e_bx);
  emulatorResults->SetBranchAddress("n"+nPrefix,&e_N);
  unpackerResults->SetBranchAddress(prefix + "_hwPt",&u_hwPt);
  unpackerResults->SetBranchAddress(prefix + "_hwEta",&u_hwEta);
  unpackerResults->SetBranchAddress(prefix + "_hwPhi",&u_hwPhi);
  unpackerResults->SetBranchAddress(prefix + "_bx",&u_bx);
  unpackerResults->SetBranchAddress("n"+nPrefix,&u_N);

  long misses = 0;
  long entries = emulatorResults->GetEntries();
  for(long i = 0; i < entries; i++)
  {
    emulatorResults->GetEntry(i);
    unpackerResults->GetEntry(i);

    int e_offset = 0;
    int u_offset = 0;

    // e_offset will always be 0, but make sure that we are not one of the
    // 1/100 events with multiple RCT BX.
    // there are 4 taus per event, 8 jets, 8 egammas, and 1 centrality
    if(e_N == 4)
    {
      e_offset = 0;
    }
    else if (e_N == 8)
    {
      e_offset = 0;
    }
    else if (e_N == 1)
    {
      e_offset = 0;
    }
    else
    {
      // there are multiple BX here and I don't know the alignment
      // skip the event
      continue;
    }

    // there are 20 taus per events, 40 jets, 40 egammas, and 5 centralities because of extra BX
    if (u_N == 5)
    {
      u_offset = 2;
    }
    else if (u_N == 20)
    {
      u_offset = 8;
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  Tin->SetBranchAddress("Padc21", &Padc[21]);
  Tin->SetBranchAddress("Padc22", &Padc[22]);
  Tin->SetBranchAddress("Padc23", &Padc[23]);
  Tin->SetBranchAddress("Padc24", &Padc[24]);
  Tin->SetBranchAddress("Padc25", &Padc[25]);
  Tin->SetBranchAddress("Padc26", &Padc[26]);
  Tin->SetBranchAddress("Padc27", &Padc[27]);
  Tin->SetBranchAddress("Padc28", &Padc[28]);
  Tin->SetBranchAddress("Padc29", &Padc[29]);
  Tin->SetBranchAddress("Padc30", &Padc[30]);
  Tin->SetBranchAddress("Padc31", &Padc[31]);

  Tin->SetBranchAddress("S83028", &S83028);
  Tin->SetBranchAddress("S8200", &S8200);
  Tin->SetBranchAddress("Clk0", &Clk0);
  Tin->SetBranchAddress("Clk1", &Clk1);
  Tin->SetBranchAddress("Clk2", &Clk2);
  Tin->SetBranchAddress("Clk3", &Clk3);

  Tin->SetBranchAddress("Pdc38",  &Pdc38);
  Tin->SetBranchAddress("Pdc39",  &Pdc39);
  Tin->SetBranchAddress("Pdc310", &Pdc310);
  Tin->SetBranchAddress("Pdc311", &Pdc311);

  Tin->SetBranchAddress("Qadc9",  &Qadc9);
  Tin->SetBranchAddress("Tdc019", &Tdc019);
  Tin->SetBranchAddress("Pdc313", &Pdc313);
  Tin->SetBranchAddress("Pdc315", &Pdc315);
  Tin->SetBranchAddress("Qadc8",  &Qadc8);
  Tin->SetBranchAddress("Tdc018", &Tdc018);
  Tin->SetBranchAddress("Qadc10", &Qadc10);
  Tin->SetBranchAddress("Tdc020", &Tdc020);

  int nEvents = Tin->GetEntries();
  cout << "... Processing nEvents = " << nEvents << endl;

  // Loop over events
  for (int i=0; i<nEvents; i++) {
    Tin->GetEvent(i);
    int iSis00 = (int) Sis00;

    // Calculate pedestal-subtracted PMT QADC values
    for (int j=0; j<8; j++) {
      pmt[j] = ((double) Qadc[j]) - pedQadc[j];
    }

    // Calculate pedestal-subtracted MWPC cathode PADC values
    for (int j=0; j<32; j++) {
      cathodeEast[j] = ((double) Pdc2[j]) - pedPdc2[j];
      cathodeWest[j] = ((double) Padc[j]) - pedPadc[j];
    }

    // Calculate pedestal-subtracted MWPC Anode PADC values
    //Took this out for now. We are comparing against a cut instead...
    //AnodeE = ((double) Pdc30) - pedPdc30;
    //AnodeW = ((double) Pdc34) - pedPdc34;

    // UCN monitor events
    bool UCNMonitorTrigger = false;
    if ( (iSis00 == 260) || (iSis00 == 516) || (iSis00 == 1028) || (iSis00 == 2052) ) {
      UCNMonitorTrigger = true;
    }

    // Events with a muon hit
    bool muonHitEast = false;
    bool muonHitWest = false;

//.........这里部分代码省略.........
   TH1D *muEta = new TH1D("muEta","muEta", 100,-3.0,3.0);
   muEta->Sumw2();


   TH1D *muPtWT = new TH1D("muPtWT","muPtWT", 100,0,10.0);
   muPtWT->Sumw2();
   
   TH1D *muEtaWT = new TH1D("muEtaWT","muEtaWT", 100,-3.0,3.0);
   muEtaWT->Sumw2();
   

   TH1D *muPhi = new TH1D("muPhi","muPhi", 100,-5.0,5.0);
   muPhi->Sumw2();

   TH1D *muDxy = new TH1D("muDxy","muDxy", 100,-30.0,30.0);
   TH1D *muDz = new TH1D("muDz","muDz", 100,-50.0,50.0);



   TH1D *VtxX = new TH1D("VtxX","VtxX", 100,-30.0,30.0);
   TH1D *VtxY = new TH1D("VtxY","VtxY", 100,-30.0,30.0);
   TH1D *VtxZ = new TH1D("VtxZ","VtxZ", 100,-50.0,50.0);


   TH1D *genmuPt = new TH1D("genmuPt","genmuPt", 100,0,10.0);
   genmuPt->Sumw2();
   TH1D *genmuEta = new TH1D("genmuEta","genmuEta", 100,-3.0,3.0);
   genmuEta->Sumw2();
   TH1D *genmuPhi = new TH1D("genmuPhi","genmuPhi", 100,-5.0,5.0);
   genmuPhi->Sumw2();

   

   Long64_t nentries = SingleWTPATMuonTree->GetEntries();
   
   for (Long64_t i=0; i<nentries;i++) {
     
     SingleWTPATMuonTree->GetEntry(i);
     
     if(i%25000==0){cout<<i<<endl;}

     VtxX->Fill(PAT_vertexXEv); 
     VtxY->Fill(PAT_vertexYEv); 
     VtxZ->Fill(PAT_vertexZEv);

     
     if(PAT_GenParSize>0){
              
       for (Int_t genpar1=0; genpar1<PAT_GenParSize;genpar1++){
	 
	 int genmuIn=0;

	 if(TMath::Abs(PAT_GenParId[genpar1])==13 && IsAccept(PAT_GenParPt[genpar1], PAT_GenParEta[genpar1]) ){genmuIn=1;}else{genmuIn=0;}
	
	   if(TMath::Abs(PAT_GenParId[genpar1])==13 && genmuIn==1 ){
	 genmuPt->Fill(PAT_GenParPt[genpar1]);
	 genmuEta->Fill(PAT_GenParEta[genpar1]);
	 genmuPhi->Fill(PAT_GenParPhi[genpar1]);
	 }
       }
     
     }//gen par size >0


     if(WTPATMuSize>0){
       //cout<<" PAT mu size "<<WTPATMuSize<<endl;

//------------------------------------------------------------------------------
// GetBoxPopulation
//------------------------------------------------------------------------------
void GetBoxPopulation(TString sample,
		      float*  box1,
		      float*  box2,
		      float&  region1_box1_yield,
		      float&  region1_box2_yield,
		      float&  region2_box1_yield,
		      float&  region2_box2_yield)
{
  TFile* file = new TFile(_inputdir + sample + ".root", "read");

  TTree* tree = (TTree*)file->Get("latino");

  float eventW;
  float mva;
  float nbjet;
  float njet;

  tree->SetBranchAddress("eventW",         &eventW);
  tree->SetBranchAddress("mva_" + _signal, &mva);
  tree->SetBranchAddress("nbjet20loose",   &nbjet);
  tree->SetBranchAddress("njet",           &njet);

  region1_box1_yield = 0;
  region1_box2_yield = 0;
  region2_box1_yield = 0;
  region2_box2_yield = 0;

  Long64_t nentries = tree->GetEntries();

  for (Long64_t ievt=0; ievt<nentries; ievt++) {

    tree->GetEntry(ievt);

    bool reject_region1 = (mva < _cut-0.2 || mva > _cut);
    bool reject_region2 = (mva < _cut-0.4 || mva > _cut-0.2);

    bool reject_box1 = ((box1[njmin] > -1 && njet  < box1[njmin]) ||
			(box1[njmax] > -1 && njet  > box1[njmax]) ||
			(box1[nbmin] > -1 && nbjet < box1[nbmin]) ||
			(box1[nbmax] > -1 && nbjet > box1[nbmax]));

    bool reject_box2 = ((box2[njmin] > -1 && njet  < box2[njmin]) ||
			(box2[njmax] > -1 && njet  > box2[njmax]) ||
			(box2[nbmin] > -1 && nbjet < box2[nbmin]) ||
			(box2[nbmax] > -1 && nbjet > box2[nbmax]));

    if (sample != "01_Data"){

	if (!reject_region1) region1_box1_yield += eventW;
        if (!reject_region2) region2_box1_yield += eventW;

    }

    else{

    if (!reject_region1 && !reject_box1) region1_box1_yield += eventW;
    if (!reject_region1 && !reject_box2) region1_box2_yield += eventW;
    if (!reject_region2 && !reject_box1) region2_box1_yield += eventW;
    if (!reject_region2 && !reject_box2) region2_box2_yield += eventW;

    }

  }

  if(sample != "01_Data") cout << sample << " -- " << region1_box1_yield << " -- " << region2_box1_yield << endl;
}

//------------------------------------------------------------------------------
// GetPdfQcdSyst
//------------------------------------------------------------------------------
void GetPdfQcdSyst(TString sample)
{

  _thunc[0] = 0; 
  _thunc[1] = 0;
  _thunc[2] = 0;
  _thunc[3] = 0;
  _thunc[4] = 0;  
  _thunc[5] = 0;


  TFile* file = new TFile(_inputdir + sample + ".root", "read");

  TTree* tree = (TTree*)file->Get("latino");

  float          mva;
  float          njet;
  vector<float> *LHEweight = 0;

  tree->SetBranchAddress("mva_" + _signal, &mva);
  tree->SetBranchAddress("njet",           &njet);
  tree->SetBranchAddress("LHEweight",      &LHEweight );

  TH1D* h_pdfsum_gen = (TH1D*)file->Get("h_pdfsum");
  TH1D* h_qcdsum_gen = (TH1D*)file->Get("h_qcdsum");

  TH1D* h_pdfsum_rec = (TH1D*)h_pdfsum_gen->Clone("h_pdfsum_rec");
  TH1D* h_qcdsum_rec = (TH1D*)h_pdfsum_gen->Clone("h_qcdsum_rec");

  h_pdfsum_rec->Reset();
  h_qcdsum_rec->Reset();

  int nbinpdf = h_pdfsum_gen->GetNbinsX();
  int nbinqcd = h_qcdsum_gen->GetNbinsX();


  // Loop over the tree
  //----------------------------------------------------------------------------
  Long64_t nentries = tree->GetEntries();

  for (Long64_t i=0; i<nentries; i++) {

    tree->GetEntry(i);

    if (mva < _cut) continue;

    if (njet < 1) continue;

    for (int a=0; a<nbinpdf; a++) h_pdfsum_rec->Fill(a, LHEweight->at(a+9));
    for (int a=0; a<nbinqcd; a++) h_qcdsum_rec->Fill(a, LHEweight->at(a));
  }


  // Produce the QCD uncertainties
  //----------------------------------------------------------------------------
  float qcdratio_gen_up   = h_qcdsum_gen->GetBinContent(9) / h_qcdsum_gen->GetBinContent(1);
  float qcdratio_gen_down = h_qcdsum_gen->GetBinContent(5) / h_qcdsum_gen->GetBinContent(1);

  float qcdratio_rec_up   = h_qcdsum_rec->GetBinContent(9) / h_qcdsum_rec->GetBinContent(1);
  float qcdratio_rec_down = h_qcdsum_rec->GetBinContent(5) / h_qcdsum_rec->GetBinContent(1);


  // Produce the PDF uncertainties
  //----------------------------------------------------------------------------
  float pdf_gen_mean = h_pdfsum_gen->Integral() / nbinpdf;
  float pdf_rec_mean = h_pdfsum_rec->Integral() / nbinpdf;

  float pdf_gen_mean_sq = 0;
  float pdf_rec_mean_sq = 0;

  for (int a=1; a<=nbinpdf; a++)
    {
      pdf_gen_mean_sq += (pow(h_pdfsum_gen->GetBinContent(a), 2) / nbinpdf);
      pdf_rec_mean_sq += (pow(h_pdfsum_rec->GetBinContent(a), 2) / nbinpdf);
    }

  float pdf_gen_sd = sqrt(pdf_gen_mean_sq - pow(pdf_gen_mean, 2)); 
  float pdf_rec_sd = sqrt(pdf_rec_mean_sq - pow(pdf_rec_mean, 2)); 

  float pdf_gen_ratio = 1 + (pdf_gen_sd / pdf_gen_mean);
  float pdf_rec_ratio = 1 + (pdf_rec_sd / pdf_rec_mean);


  // Print the uncertainties
  //----------------------------------------------------------------------------
  //printf("\n %s\n", sample.Data());
  //printf("---------------------------------------\n");
  //printf(" QCD up   -- xs = %5.3f -- acc = %5.3f\n", qcdratio_gen_up,   qcdratio_rec_up   / qcdratio_gen_up);
  //printf(" QCD down -- xs = %5.3f -- acc = %5.3f\n", qcdratio_gen_down, qcdratio_rec_down / qcdratio_gen_down);
  //printf(" PDF      -- xs = %5.3f -- acc = %5.3f\n", pdf_gen_ratio,     pdf_rec_ratio     / pdf_gen_ratio);

  _thunc[0] = qcdratio_gen_up                    ; 
  _thunc[1] = qcdratio_rec_up/qcdratio_gen_up    ;
  _thunc[2] = qcdratio_gen_down                  ;
  _thunc[3] = qcdratio_rec_down/qcdratio_gen_down;
  _thunc[4] = pdf_gen_ratio                      ;  
  _thunc[5] = pdf_rec_ratio/pdf_gen_ratio        ; 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	//CUTS__________________________
	double lowLimitPt    = ptBinStart[0];
	double hiLimitPt     = ptBinEnd[4];
	double lowLimitEta   = etaBinStart[0];
	double hiLimitEta    = etaBinEnd[3];
	double lowLimitMass  = massBinStart[0];
	double hiLimitMass   = massBinEnd[4];
	
	//hiLimitEta = 0;
	
	if (fullEta)
	{
		assert(hiLimitEta > 0);
		assert(lowLimitEta < 0);
	}
	
	
	cout << "Pt  between " << lowLimitPt << " and " << hiLimitPt << endl;
	cout << "M   between " << lowLimitMass << " and " << hiLimitMass << endl;
	cout << "Eta between " << lowLimitEta << " and " << hiLimitEta << endl;
	//______________________________
	
	//*
	
	//   ======================================================================
	//============================================================================
	//START ANALYSIS==============================================================
	//============================================================================
	//   ======================================================================
	//*
	cout << "\n";
	cout << "<----STARTING ANALYSIS---->" << endl;
  cout << "\n";
	cout << pairTree->GetEntries() << " pairs to analyze" << endl;
	
	
	
	int blueFillNo;
	int yellowFillNo;
	
	int phiSRbin;
	
	TLorentzVector sum;
	TLorentzVector sumY;
	TLorentzVector sumB;
	
	//random number for randomizing spin.
	//set seed to zero to gaurenty unique numbers each time.
	TRandom3 r;
	r.SetSeed(0);
	
	int rand5 = 0;
	int rand6 = 0;
	int rand9 = 0;
	int rand10 = 0;
	
	
	int totalPairsFinal = 0;
	int blueD = 0;
	int blueU = 0;
	int yellowD = 0;
	int yellowU = 0;
	
	int pionStarNumber = 0;
	int test = 0;
	

int main (int argc, char** argv)
{ 
  std::string inputFolder = argv[1];
  std::string outputFile = argv[2];
  std::string leptonName = argv[3];
  std::string inputFile = argv[4];

  std::cout<<"file: "<<(inputFile).c_str()<<std::endl;
  TFile * fS = new TFile((inputFolder+leptonName+"/WWTree_"+inputFile+".root").c_str());
  TTree * inputTree = (TTree *)fS->Get("otree");

  int run;
  int event;
  int lumi;
  int njets;
  int nPV;
  int issignal;
  float pfMET;
  float pfMET_Phi;
  float l_pt;
  float l_eta;
  float l_phi;
  float l_e;
  float ungroomed_jet_pt;
  float ungroomed_jet_eta;
  float ungroomed_jet_phi;
  float ungroomed_jet_e;
  float jet_mass_pr;
  float jet_mass_so;
  float jet_tau2tau1;
  float v_pt;
  float v_eta;
  float v_phi;
  float v_mt;
  float mass_lvj_type0;
  int nBTagJet_medium;
  float jet2_pt;
  float jet2_btag;
  float jet3_pt;
  float jet3_btag;

  inputTree->SetBranchAddress("run", &run);
  inputTree->SetBranchAddress("event", &event);
  inputTree->SetBranchAddress("lumi", &lumi);
  inputTree->SetBranchAddress("njets", &njets);
  inputTree->SetBranchAddress("nPV", &nPV);
  inputTree->SetBranchAddress("issignal", &issignal);
  inputTree->SetBranchAddress("pfMET", &pfMET);
  inputTree->SetBranchAddress("pfMET_Phi", &pfMET_Phi);
  inputTree->SetBranchAddress("l_pt", &l_pt);
  inputTree->SetBranchAddress("l_eta", &l_eta);
  inputTree->SetBranchAddress("l_phi", &l_phi);
  inputTree->SetBranchAddress("l_e", &l_e);
  inputTree->SetBranchAddress("ungroomed_jet_pt", &ungroomed_jet_pt);
  inputTree->SetBranchAddress("ungroomed_jet_eta", &ungroomed_jet_eta);
  inputTree->SetBranchAddress("ungroomed_jet_phi", &ungroomed_jet_phi);
  inputTree->SetBranchAddress("ungroomed_jet_e", &ungroomed_jet_e);
  inputTree->SetBranchAddress("jet_mass_pr", &jet_mass_pr);
  inputTree->SetBranchAddress("jet_mass_so", &jet_mass_so);
  inputTree->SetBranchAddress("jet_tau2tau1", &jet_tau2tau1);
  inputTree->SetBranchAddress("v_pt", &v_pt);
  inputTree->SetBranchAddress("v_eta", &v_eta);
  inputTree->SetBranchAddress("v_phi", &v_phi);
  inputTree->SetBranchAddress("v_mt", &v_mt);
  inputTree->SetBranchAddress("mass_lvj_type2", &mass_lvj_type0);
  inputTree->SetBranchAddress("nBTagJet_medium", &nBTagJet_medium);
  inputTree->SetBranchAddress("jet2_pt", &jet2_pt);
  inputTree->SetBranchAddress("jet2_btag", &jet2_btag);
  inputTree->SetBranchAddress("jet3_pt", &jet3_pt);
  inputTree->SetBranchAddress("jet3_btag", &jet3_btag);

  //---------output tree----------------
  TFile* outROOT = TFile::Open((std::string("output/output_synch_")+leptonName+std::string("/")+std::string("WWTree_")+outputFile+std::string(".root")).c_str(),"recreate");
  outROOT->cd();
  TTree* outTree = new TTree("otree", "otree");
  outTree->SetDirectory(0);
  setOutputTreeSynch *WWTree = new setOutputTreeSynch(outTree);
  
  //---------start loop on events------------
  for (Long64_t jentry=0; jentry<inputTree->GetEntries();jentry++) {

    inputTree->GetEntry(jentry);

    WWTree->initializeVariables(); //initialize all variables

    if(jentry % 1000 == 0)    
      cout << "read entry: " << jentry << endl;

    WWTree->issignal = issignal;

    //save event variables
    WWTree->run   = run;
    WWTree->event = event;
    WWTree->lumi = lumi;
    WWTree->njets = njets;
    WWTree->nPV  = nPV;
    
    WWTree->l_pt  = l_pt;
    WWTree->l_eta = l_eta;
    WWTree->l_phi = l_phi;	
//.........这里部分代码省略.........

int main (int argc, char **argv) {

  TFile *tf = TFile::Open("root/AnalysisOut.root");
  TTree *tree = (TTree*)tf->Get("AnalysisTree");

  RooWorkspace *w = new RooWorkspace("w","w");

  // tree variables
  Double_t  t_mass;
  Double_t  t_KST1_m;
  Double_t  t_KST2_m;
  int       t_itype;
  ULong64_t t_eventNumber;
  Bool_t    t_pass_bdt;
  Bool_t    t_pass_pid;
  TString   *t_evname = 0;

  // tree branches
  tree->SetBranchAddress("B_s0_DTF_B_s0_M",      &t_mass        );
  tree->SetBranchAddress("B_s0_DTF_KST1_M",      &t_KST1_m   );
  tree->SetBranchAddress("B_s0_DTF_KST2_M",      &t_KST2_m   );
  tree->SetBranchAddress("itype",                &t_itype       );
  tree->SetBranchAddress("eventNumber",          &t_eventNumber );
  tree->SetBranchAddress("pass_bdt",             &t_pass_bdt    );
  tree->SetBranchAddress("pass_pid",             &t_pass_pid    );
  tree->SetBranchAddress("evname",               &t_evname      );

  // observables
  RooRealVar *mass        = new RooRealVar("B_s0_DTF_B_s0_M", "m(K#piK#pi)", 5000., 5800.);
  RooRealVar *KST1_mass   = new RooRealVar("B_s0_DTF_KST1_M", "m(K#pi)"    , 750, 1700 );
  RooRealVar *KST2_mass   = new RooRealVar("B_s0_DTF_KST2_M", "m(K#pi)"    , 750, 1700 );
  RooRealVar *itype       = new RooRealVar("itype","itype",-10,10); itype->removeRange();
  RooRealVar *eventNumber = new RooRealVar("eventNumber","eventNumber",0,10e10);


  // map to store datasets
  map<int, RooDataSet*> dsets;

  for ( int ev=0; ev<tree->GetEntries(); ev++ ) {
    tree->GetEntry(ev);

    if ( t_itype >= 0 && !t_pass_bdt ) continue;
    if ( t_itype >= 0 && !t_pass_pid ) continue;

    // if entry not in map then add it
    if ( dsets.find( t_itype ) == dsets.end() ) {
      dsets[t_itype] = new RooDataSet( *t_evname, *t_evname, RooArgSet(*mass,*KST1_mass,*KST2_mass,*itype,*eventNumber) );
    }

    // add to dataset
    mass->setVal( t_mass );
    KST1_mass->setVal( t_KST1_m );
    KST2_mass->setVal( t_KST2_m );
    itype->setVal( t_itype );
    eventNumber->setVal( t_eventNumber );
    dsets[t_itype]->add( RooArgSet(*mass,*KST1_mass,*KST2_mass,*itype,*eventNumber) );

    if ( ev%10000==0 ) {
      cout << ev << "/" << tree->GetEntries() << endl;
      cout << "\t" << t_itype << " " << *t_evname << " " << t_eventNumber << " " << t_mass << endl;
    }
  }

  for ( map<int, RooDataSet*>::iterator it = dsets.begin(); it != dsets.end(); it++ ) {

    cout << "\t" << it->first << " " << it->second->GetName() << " " << it->second->numEntries() << endl;
    w->import(*it->second);

  }

  tf->Close();

  delete t_evname;

  w->writeToFile("tmp/DataSets.root");
  return 0;
}