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

plot_warnmap2D( string warnmapfile , bool writeplots = true )
{
  gStyle->SetOptStat(0);

  /* Read warnmap as tree */
  TTree *twarn = new TTree();
  twarn->ReadFile( warnmapfile.c_str(), "sector/I:y:z:status" );

  // Declaration of leaf types
  Int_t           sector;
  Int_t           y;
  Int_t           z;
  Int_t           status;

  // List of branches
  TBranch        *b_sector;   //!
  TBranch        *b_y;   //!
  TBranch        *b_z;   //!
  TBranch        *b_status;   //!

  twarn->SetBranchAddress("sector", &sector, &b_sector);
  twarn->SetBranchAddress("y", &y, &b_y);
  twarn->SetBranchAddress("z", &z, &b_z);
  twarn->SetBranchAddress("status", &status, &b_status);



  /* Create 2D histograms to visualize warnmap */
  TH2I* h_warnmap[8];
  int zbins = 0;
  int ybins = 0;

  /* Create histogram for each sector */
  for( int sector = 0; sector < 8; sector++ )
    {
      if(sector <6) // PbSc sector
	{
	  zbins = 72;
	  ybins = 36;
	}
      else // PbGl sector
	{
	  zbins = 96;
	  ybins = 48;
	}

      TString name("h_warnmap_sector_");
      name += sector;

      TString title("Warnmap sector ");
      title += sector;

      cout << name <<endl;
      cout << title <<endl;

      h_warnmap[sector] = new TH2I( name, title, zbins, 0, zbins, ybins, 0, ybins );
      h_warnmap[sector]->GetXaxis()->SetTitle("z [cm]");
      h_warnmap[sector]->GetYaxis()->SetTitle("y [cm]");

    } //sector

  /* Loop over warnmap and fill status of towers into map */
  Long64_t nentries = twarn->GetEntriesFast();

  for ( Long64_t jentry = 0; jentry < nentries; jentry++ )
    {
      twarn->GetEntry(jentry);

      //cout << sector << " " << y << " " << z << " " << status << endl;

      // if status > 50 (e.g. 'dead' or 'uncalibrated', set it to 39 for nicer plotting
      if ( status > 50 )
	status = 39;

      h_warnmap[sector]->Fill( z , y , status );

    }

  /* build base filename for plots */
  std::size_t pos = warnmapfile.find("/");
  string filename_cut = warnmapfile.substr( pos+1 );
  (filename_cut.erase( filename_cut.length()-4 ,4 ));//.erase(0,25);


  /* plot warnmaps */
  for( int sector = 0; sector < 8; sector++ )
    {
      TCanvas *c1 = new TCanvas();
      (h_warnmap[sector])->Draw("colz");

      if ( writeplots )
	{
	  TString filename("plots-warnmap-2D/warnmap2D_");
	  filename+=filename_cut;
	  filename+="_sector_";
	  filename+=sector;
	  filename+=".eps";

	  TString filenamep("plots-warnmap-2D/warnmap2D_");
	  filenamep+=filename_cut;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	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 runsProcessed = 0;
	
	double currentFillNo = 0;
	double currentRunNo = 0;

	
	for (int iPair = pionStarNumber; iPair < pairTree->GetEntries(); iPair++)
	{
		
		if (iPair%10000 == 0) {cout << "processing pair number " << iPair << endl;}
		//if (iPair == pionStarNumber+300000){break;}
		
		pairTree->GetEntry(iPair);
		
		//if (runsProcessed > 4){break;} 
		
		
		if (pair1->runInfo().beamFillNumber(1) != currentFillNo && currentFillNo != 0 && currentFillNo != 16427)
		//if (pair1->runInfo().runId() != currentRunNo && currentRunNo != 0)
		{
			runsProcessed++;

			cout << "runsProcessed " << runsProcessed << endl;
			
			double* asym = new double();
			double* asymE = new double();
			calcAsyms(hAut, hNumberUp, hNumberDown, polOfBin, asym, asymE);
						
			fillAsyms.push_back(*asym);
			fillAsymsE.push_back(*asymE);
			fillVec.push_back(runsProcessed);
			
			hNumberUp->Reset();
			hNumberDown->Reset();
			hAut->Reset();
			
			for (int i=0; i<binNumber; i++)
			{
				polOfBin[i]->Reset();
			}
			
			
		}
		else if (currentRunNo == 16427)
		{

//.........这里部分代码省略.........
            reader->BookMVA( methodName, weightfile );
        }
    }
    
    // Prepare input tree (this must be replaced by your data source)
    // in this example, there is a toy tree with signal and one with background events
    // we'll later on use only the "signal" events for the test in this example.
    //
    
    /*//TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/data_merged_ccqe_addpid_ver3noveract_pid1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pm_merged_ccqe_tot_addpid_ver3noveract_pid1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pmbar_merged_ccqe_addpid_ver3noveract_pid1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/ingrid_merged_nd3_ccqe_tot_addpid_ver3noveract_pid1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/wall_merged_ccqe_tot_addpid_ver3noveract_pid1pres.root";
    //check genie signal
    TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/genie_merged_ccqe_coh_addpid_ver3noveract_pid1pres.root";*/
    
    //add for event with more than 2 track
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/data_merged_ccqe_processl2trk.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pm_merged_ccqe_tot_processl2trk.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pmbar_merged_ccqe_processl2trk.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/ingrid_merged_nd3_ccqe_tot_processl2trk.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/wall_merged_ccqe_tot_processl2trk.root";
    //check genie signal
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/genie_merged_ccqe_coh_processl2trk.root";
    
    //for correct cc1pres pid
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/data_merged_ccqe_correct1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pm_merged_ccqe_tot_correct1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pmbar_merged_ccqe_correct1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/ingrid_merged_nd3_ccqe_tot_correct1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/wall_merged_ccqe_tot_correct1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/genie_merged_ccqe_coh_correct1pres.root";
    //for correct prange/pang_t and additional information
    //TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/data_merged_ccqe_addpidFF.root";
   // TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/pm_merged_ccqe_tot_addpidFF.root";
    //TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/pmbar_merged_ccqe_addpidFF.root";
    //TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/ingrid_merged_nd3_ccqe_tot_addpidFF.root";
    //TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/wall_merged_ccqe_tot_addpidFF.root";
    //TString fname = "/home/cvson/cc1picoh/dataProcess/fix20150420/genie_merged_ccqe_tot_addpidFF.root";
    
    //
    std::cout << "--- Selecting data sample" << std::endl;
    TFile *pfile = new TFile(fname,"update");
    TTree* theTree = (TTree*)pfile->Get("tree");
    theTree->SetBranchAddress( "mumucl", &mumucl );
    theTree->SetBranchAddress( "pmucl", &pmucl );
    theTree->SetBranchAddress( "pang_t", &pang_t );
    theTree->SetBranchAddress( "muang_t", &muang_t );
    theTree->SetBranchAddress( "veract", &veract );
    theTree->SetBranchAddress( "ppe", &ppe);
    theTree->SetBranchAddress( "mupe", &mupe);
    theTree->SetBranchAddress( "range", &range);
    theTree->SetBranchAddress( "coplanarity", &coplanarity);
    theTree->SetBranchAddress( "opening", &opening);
    
 
    
    Int_t Ntrack;
    theTree->SetBranchAddress( "Ntrack", &Ntrack );
    
    
    Float_t pid1pres;
    TBranch *bpid1pres = theTree->Branch("pid1pres",&pid1pres,"pid1pres/F");
    
    
    
    std::vector<Float_t> vecVar(9); // vector for EvaluateMVA tests
    
    Long64_t nentries = theTree->GetEntriesFast();
    Long64_t iprintProcess = Long64_t(nentries/100.);
    
    std::cout << "--- Processing: " << nentries << " events" << std::endl;
    
    TStopwatch sw;
    sw.Start();
    for (Long64_t ievt=0; ievt<nentries;ievt++) {
        
        if (ievt%iprintProcess == 0) cout<<"Processing "<<int(ievt*100./nentries)<<"% of events"<<endl;
        theTree->GetEntry(ievt);
        Float_t pid_tem;
        if (Use["BDT"])      {
            //if (Ntrack!=2) pid_tem = -999;//change here
            if (Ntrack<2) pid_tem = -999;
            else pid_tem = reader->EvaluateMVA("BDT method");
        }
        pid1pres = pid_tem;
        bpid1pres->Fill();
        
    }
    theTree->Write();
    delete pfile;
    // Get elapsed time
    sw.Stop();
    std::cout << "--- End of event loop: "; sw.Print();
    
    delete reader;
    
    std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
} 

//.........这里部分代码省略.........
    tree->SetBranchAddress("muNeg_Trigger2", &muNeg_Trigger2);
    tree->SetBranchAddress("muNeg_Trigger21", &muNeg_Trigger21);
    tree->SetBranchAddress("muNeg_Trigger22", &muNeg_Trigger22);
    //====================================Gen Variables=========================================================
    //Gen Jpsi Variables
    gentree->SetBranchAddress("GenJpsiMass",   &GenJpsiMass);
    gentree->SetBranchAddress("GenJpsiPt",     &GenJpsiPt);
    gentree->SetBranchAddress("GenJpsiPhi",    &GenJpsiPhi);
    gentree->SetBranchAddress("GenJpsiRap",    &GenJpsiRap);
    gentree->SetBranchAddress("GenJpsiEta",    &GenJpsiEta);
    gentree->SetBranchAddress("GenJpsiPx",     &GenJpsiPx);
    gentree->SetBranchAddress("GenJpsiPy",     &GenJpsiPy);
    gentree->SetBranchAddress("GenJpsiPz",     &GenJpsiPz);
    //muon variable
    gentree->SetBranchAddress("GenmuPosPx",    &GenmuPosPx);
    gentree->SetBranchAddress("GenmuPosPy",    &GenmuPosPy);
    gentree->SetBranchAddress("GenmuPosPz",    &GenmuPosPz);
    gentree->SetBranchAddress("GenmuPosEta",   &GenmuPosEta);
    gentree->SetBranchAddress("GenmuPosPhi",   &GenmuPosPhi);
    gentree->SetBranchAddress("GenmuNegPx",    &GenmuNegPx);
    gentree->SetBranchAddress("GenmuNegPy",    &GenmuNegPy);
    gentree->SetBranchAddress("GenmuNegPz",    &GenmuNegPz);
    gentree->SetBranchAddress("GenmuNegEta",   &GenmuNegEta);
    gentree->SetBranchAddress("GenmuNegPhi",   &GenmuNegPhi);
    //gentree->SetBranchAddress("GenvertexZ",&GenvertexZ);

    //====================================================== Gen tree loop ================================================
    int NAccep=0;
    int nGenEntries=gentree->GetEntries();
    cout<<" Total Entries in GenLevel Tree for pT range: "<<fileName[0]<<"  "<<   nGenEntries<< " ==============="<<endl;

    for(int i=0; i< nGenEntries; i++)  {        
      //cout<<"i : "<<i<<endl;
      if(!(gentree->GetEntry(i))) continue;
      //cout<<" gentree ("<<i<<")"<<endl;
      //Only printing 
      if(i%1000000==0){
        cout<<" processing record "<<i<<"/"<<nGenEntries<<endl;
        //cout<<" Mass "<< GenJpsiMass<< " pT "<< GenJpsiPt << " Y " <<GenJpsiRap<<endl;
      }

      bool GenPosIn=0, GenNegIn=0;
      GenmuPosPt= TMath::Sqrt(GenmuPosPx*GenmuPosPx + GenmuPosPy*GenmuPosPy); 
      GenmuNegPt= TMath::Sqrt(GenmuNegPx*GenmuNegPx + GenmuNegPy*GenmuNegPy); 

      if(IsAccept(GenmuPosPt, GenmuPosEta)) {GenPosIn=1;}
      if(IsAccept(GenmuNegPt, GenmuNegEta)) {GenNegIn=1;}

      int AccJpsi = 0;
      //if(GenJpsiPt < 6.5) continue;
      if((GenJpsiPt >= minPt && GenJpsiPt <= maxPt && fabs(GenJpsiRap) >= minRap && TMath::Abs(GenJpsiRap) <= maxRap && GenPosIn == 1 && GenNegIn == 1)
          && GenJpsiMass > 2.0 && GenJpsiMass < 4.0
          ) {AccJpsi = 1;}

      if((GenPosIn ==1 && GenNegIn==1)) NAccep++;

      //cout<<"1. GenJpsiPt : "<<GenJpsiPt<<", GenJpsiEta : "<<GenJpsiEta<<", GenJpsiRap : "<<GenJpsiRap<<", |GenJpsiPsi| : "<<TMath::Abs(GenJpsiPsi)<<endl;
      //cout<<"1. GenPosIn : "<<GenPosIn<<", GenNegIn : "<<GenNegIn<<endl;
      double vars = 0.0, bin1 = 0.0, bin2 = 0.0;
      if(iSpec == 0) vars = GenJpsiPt;
      if(iSpec == 1) vars = fabs(GenJpsiRap);
      for(int j = 0; j < nBins; j++){
        if(iSpec == 0){
          bin1 = pt_bound[j]; bin2 = pt_bound[j+1];
          if( (AccJpsi==1) && (vars >= bin1 && vars < bin2) && fabs(GenJpsiRap) >= minRap && TMath::Abs(GenJpsiRap) <= maxRap && TMath::Abs(GenJpsiRap) >= minRap) {
            hGenDiMuon[j]->Fill(GenJpsiMass);

int main(int argc, char* argv[])
{
  
    printf("Before tool\n");
    BCHTool::BCHCleaningToolRoot thebchtool("thebchtool");
    thebchtool.SetResponseVsBCHFile("BCHCleaningTool/share/BCH_Response.root");

    ///TileTripReader initialisation and tests outside BCH tool
    
    //prepare a TileTripReader instance
    Root::TTileTripReader *ttr=new Root::TTileTripReader("thetripreader");
    ttr->setTripFile("TileTripReader/data/CompleteTripList_2011-2012.root");
    ttr->setVerbosity(0);
    
    //here we ask the TTR directly to give us the answer of a jet that is on the edge
    //the edges haven't yet been changed, so the answer should be no.
    bool isInDeadRegionBeforeEdgeChange = ttr->checkEtaPhi(211787,25000,0.1,0.35);
    cout << "In dead region, before edge change? " << isInDeadRegionBeforeEdgeChange << endl;

    //Set the edges
    ttr->m_LBOffsets.eta1=0.1;
    ttr->m_LBOffsets.eta2=0.1;
    ttr->m_LBOffsets.phi1=TMath::Pi()/32.;
    ttr->m_LBOffsets.phi2=TMath::Pi()/32.;
    
    //here we ask the TTR directly to give us the answer of a jet that is on the edge
    //the edges have now been changed, so the answer should be yes.
    bool isInDeadRegionAfterEdgeChange = ttr->checkEtaPhi(211787,25000,0.1,0.35);
    //NOTE: THIS IS WHERE IT DOES NOT GIVE ME THE RIGHT ANSWER!!! it should be TRUE below
    cout << "In dead region, after edge change? " << isInDeadRegionAfterEdgeChange << endl;
    
    ///PRW initialisation and tests
    
    //prepare a PRW instance
    Root::TPileupReweighting *prw=new Root::TPileupReweighting("thepileupreweighting");
    //initialise it with a MC file, a data file and tell it what to do with unmatched data
    prw->AddConfigFile("RootFilesForTests/ExcitedQ_2000.prw.root");
    prw->AddLumiCalcFile("RootFilesForTests/ilumicalc_histograms_None_200841-204158.root"); 
    prw->SetUnrepresentedDataAction(2);
    prw->Initialize();
    
    //initialize the BCH tool (for MC at the moment)
    //this will change the edges for the offsets of the TTR
    thebchtool.InitializeTool(false, ttr, prw);
    
    //here we ask the TTR directly to give us the answer of a jet that is on the edge
    //the edges haven't yet been changed, so the answer should be no.
    isInDeadRegionAfterEdgeChange = ttr->checkEtaPhi(211787,25000,0.1,0.35);
    cout << "In dead region, after edge change from BCH tool? " << isInDeadRegionAfterEdgeChange << endl;

    ///Stress test 
    TFile* stressfulFile = TFile::Open("RootFilesForTests/stress_test.root");
    TTree* stressfulTree = (TTree*) stressfulFile->Get("physics");
   
    ///Some random q* MC 
    //TODO add a different sample here
    
    std::vector<float> * v_jets_pt=0;
    std::vector<float> * v_jets_eta=0;
    std::vector<float> * v_jets_phi=0;
    stressfulTree->SetBranchAddress("jet_AntiKt4LCTopo_pt", &v_jets_pt);
    stressfulTree->SetBranchAddress("jet_AntiKt4LCTopo_eta", &v_jets_eta);
    stressfulTree->SetBranchAddress("jet_AntiKt4LCTopo_phi", &v_jets_phi);

    for (Int_t i=0;i<stressfulTree->GetEntries();i++) {
      stressfulTree->GetEntry(i);
      for (UInt_t ijet=0; ijet<v_jets_pt->size(); ijet++) {
        
        cout <<  "Jet pt:" << v_jets_pt->at(ijet) << endl;
        cout <<  "Jet eta:" << v_jets_eta->at(ijet) << endl;
        cout <<  "Jet phi:" << v_jets_phi->at(ijet) << endl;

        ///testing the basic functionality
        
        //note: need to pass the random number from MC to the BCHCleaningTool!! One can decide at this point whether to use mu or not
        //Here I've passed the RunNumber from a MC sample
        double mu = 20;
        int randomRunNumber = prw->GetRandomRunNumber(195847, mu);
        std::cout << "Random run number: " << randomRunNumber << std::endl;
        bool isInMaskedRegionBCH = thebchtool.IsInMaskedRegion(195847, v_jets_pt->at(ijet),v_jets_eta->at(ijet),v_jets_phi->at(ijet),mu);
        cout << "In dead region from BCHCleaningTool? " << isInMaskedRegionBCH  << endl;
        //testing straight from the TTR
        ttr->setVerbosity(0);
        bool isInMaskedRegionTTR = ttr->checkEtaPhi(randomRunNumber,v_jets_pt->at(ijet),v_jets_eta->at(ijet),v_jets_phi->at(ijet));
        cout << "In dead region from TTR? " << isInMaskedRegionTTR << endl;
        
        ///testing the response estimate
        cout << "Jet pT after response correction:" << thebchtool.EstimateJetResponseNoBCH(v_jets_pt->at(ijet), 0.3) << endl;        
        
      }
    }
    
    printf("After tool\n");
    return 0;
}

void csv(TString input="tmva.csvoutput.txt", TString par1="par2", TString par2="par3", TString par3="", TString value="eventEffScaled_5") {
  std::cout << "Usage:" << std::endl
            << ".x scripts/csv.C    with default arguments" << std::endl
            << ".x scripts/csv.C(filename, par1, par2, value)" << std::endl
            << std::endl
            << "  Optional arguments:" << std::endl
            << "    filename        path to CSV file" << std::endl
            << "    par1            name of X-parameter branch" << std::endl
            << "    par2            name of Y-parameter branch (if empty, efficiency is drawn as a function of par1)" << std::endl
            << "    value           name of result (efficiency) branch" << std::endl
            << std::endl;

  TTree *tree = new TTree("data", "data");
  tree->ReadFile(input);

  gStyle->SetPalette(1);
  gStyle->SetPadRightMargin(0.14);

  TCanvas *canvas = new TCanvas("csvoutput", "CSV Output", 1200, 900);

  tree->SetMarkerStyle(kFullDotMedium);
  tree->SetMarkerColor(kRed);
  if(par2.Length() > 0) {
    //tree->Draw(Form("%s:%s", par2.Data(), par1.Data()));
    if(par3.Length() > 0)
      tree->Draw(Form("%s:%s:%s:%s", par1.Data(), par2.Data(), par3.Data(), value.Data()), "", "COLZ"); //, "", "Z");
    else
      tree->Draw(Form("%s:%s:%s", par2.Data(), par1.Data(), value.Data()), "", "COLZ"); //, "", "Z");

    TH1 *histo = tree->GetHistogram();
    if(!histo)
      return;

    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(Form("Classifier parameter %s", par2.Data()));
    if(par3.Length() > 0)
      histo->GetZaxis()->SetTitle(Form("Classifier parameter %s", par3.Data()));
    else
      histo->GetZaxis()->SetTitle("");

    if(par3.Length() == 0) {
      float x = 0;
      float y = 0;
      float val = 0;
      double maxVal = tree->GetMaximum(value);
      double minVal = tree->GetMinimum(value);

      tree->SetBranchAddress(par1, &x);
      tree->SetBranchAddress(par2, &y);
      tree->SetBranchAddress(value, &val);
      TLatex l;
      l.SetTextSize(0.03);
    
      Long64_t nentries = tree->GetEntries();
      for(Long64_t entry=0; entry < nentries; ++entry) {
        tree->GetEntry(entry);
    
        l.SetTextColor(textColor(val, maxVal, minVal));
        l.DrawLatex(x, y, Form("%.3f", val*100));
      }
    }
  }
  else {
    tree->Draw(Form("%s:%s", value.Data(), par1.Data()));
    TH1 *histo = tree->GetHistogram();
    if(!histo) 
      return;
    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(value);
  }
}

//.........这里部分代码省略.........
    hVnSub0[icent]            = new TH1D( Form("hVnSub0_c%i", icent), Form("hVnSub0_c%i", icent), NBins, 0., vnMax[norder_] );
    hVnSub0[icent]->GetXaxis()->SetTitle( Form("v_{%i}^{obs,a}", norder_) );

    hVnSub1[icent]            = new TH1D( Form("hVnSub1_c%i", icent), Form("hVnSub1_c%i", icent), NBins, 0., vnMax[norder_] );
    hVnSub1[icent]->GetXaxis()->SetTitle( Form("v_{%i}^{obs,b}", norder_) );
    
    Mult[icent]               = new TH1I( Form("Mult_c%i", icent), Form("Mult_c%i", icent), 300, 1, 1500 );
    Mult[icent]->GetXaxis()->SetTitle("Multiplicity");
	
    h2Vn2D0v1[icent]          = new TH2D( Form("h2Vn2D0v1_c%i", icent), Form("h2Vn2D0v1_c%i", icent), 2*NBins, -vnMax[norder_], vnMax[norder_], 2*NBins, -vnMax[norder_], vnMax[norder_] );
    h2Vn2D0v1[icent]->GetXaxis()->SetTitle( Form("(v_{%i,x}^{obs,a} - v_{%i,x}^{obs,b})/2", norder_, norder_) );
    h2Vn2D0v1[icent]->GetYaxis()->SetTitle( Form("(v_{%i,y}^{obs,a} - v_{%i,y}^{obs,b})/2", norder_, norder_) );
    h2Vn2D0v1[icent]->SetOption("colz");

    h2Vn2D0v1Magnitude[icent] = new TH1D( Form("h2Vn2D0v1Magnitude_c%i", icent), Form("h2Vn2D0v1Magnitude_c%i", icent), NBins, 0., vnMax[norder_] );
    h2Vn2D0v1Magnitude[icent]->GetXaxis()->SetTitle( Form("|(v_{%i,x}^{obs,a} - v_{%i,x}^{obs,b})/2|", norder_, norder_) );

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

  for(int ievent = 0; ievent < N; ievent++) {
        
    if((ievent+1)% 500000 == 0) cout<<"Processing Event "<<ievent+1<<"\t"<<100.*(ievent+1)/(double)N<<"% Completed"<<endl;
      
    tree->GetEntry(ievent);

    //-- Vertex Cut
    if(TMath::Abs(vtx) > 3.0) continue;

    //-- Calculate centbin
    if( centval > cent_max[NCENT-1]) continue;
    int icent = hCentBins.FindBin(centval)-1;

    //-- Reset raw and sumw values
    VnRaw_x_0    = 0;
    VnRaw_y_0    = 0;
    VnRaw_x_1    = 0;
    VnRaw_y_1    = 0;
    VnRaw_x_full = 0;
    VnRaw_y_full = 0;                

    sumw_0       = 0;
    sumw_1       = 0;
    sumw_full    = 0;

    evtMult_0    = 0;
    evtMult_1    = 0;
    evtMult_full = 0;

    //-- Begin analyzer histogram loops
    for(int ipt = ptBinMin; ipt <= ptBinMax; ipt++){
      for(int ieta = etaBinMin; ieta <= etaBinMax; ieta++){

	if(sumw->GetBinContent(ipt+1,ieta+1) !=0){

	  //-- Subevent 0 (eta >= 0)
	  if(etabinsDefault[ieta] >= 0){

void gastof_simple()
{
  const unsigned int run_id = 725;
  TFile* f = new TFile(Form("../../gastof_inner_run%d.root", run_id)); //FIXME
  TTree* t = (TTree*)f->Get("tdc");

  const unsigned int channels_to_probe = 32; //64;
  const double tot_min = 70., tot_max = -1.;
  const double lead_min = 6.8, lead_max = 6.95; // in \mus

  int num_hits;
  unsigned int ettt;
  double leading_edge[5000], tot[5000];
  int channel_id[5000];
  t->SetBranchStatus("*", 0);
  t->SetBranchStatus("num_measurements", 1); t->SetBranchAddress("num_measurements", &num_hits);
  t->SetBranchStatus("leading_edge", 1); t->SetBranchAddress("leading_edge", leading_edge);
  t->SetBranchStatus("channel_id", 1); t->SetBranchAddress("channel_id", channel_id);
  t->SetBranchStatus("tot", 1); t->SetBranchAddress("tot", tot);

  const unsigned int num_triggers = t->GetEntries();

  TH1D* h_tot_all_channels = new TH1D("h_tot_all_channels", "", 250, 0., 1000.);
  TH1D* h_lead_all_channels = new TH1D("h_lead_all_channels", "", 500, 0., 10.);

  double last_ettt = -1.;
  unsigned int num_overfl = 0;
  unsigned int mult[channels_to_probe];
  for (unsigned int i=0; i<num_triggers; i++) {
    t->GetEntry(i);
    for (unsigned int j=0; j<channels_to_probe; j++) {
//       lead[j].clear();    //lead[] not defined in gastof_simple
      mult[j] = 0;
    }
    if (i%10000==0) cerr << "Processing event " << i << " / " << num_triggers << endl;

    for (int j=0; j<num_hits; j++) {
      h_lead_all_channels->Fill(leading_edge[j]/1.e3, 1./num_triggers);
      if (leading_edge[j]<lead_min*1.e3 or leading_edge[j]>lead_max*1.e3) continue;
      h_tot_all_channels->Fill(tot[j]);
    }
  }

  GastofCanvas c_tot_all_chan("tot_all-channels", Form("Run %d, all channels", run_id));
  h_tot_all_channels->Draw();
  h_tot_all_channels->GetXaxis()->SetTitle("Time over threshold (ns)");
  h_tot_all_channels->GetYaxis()->SetTitle("Hits");

  {
    TPaveText* cuts = new TPaveText(0.4, 0.8, 0.5, 0.85, "ndc");
    cuts->SetTextAlign(13);
    cuts->SetTextFont(43);
    cuts->SetTextSize(18);
    cuts->SetFillColor(kWhite);
    cuts->SetLineColor(kWhite);
    //cuts->AddText(Form("ToT > %.1f ns", tot_min));
    cuts->AddText(Form("Lead. edge #in [%.1f-%.2f]  #mus", lead_min, lead_max));
    cuts->Draw("same");
  }

  c_tot_all_chan.Pad()->SetLogy();
  {
    TLine* line = new TLine(tot_min, 0., tot_min, h_tot_all_channels->GetMaximum()*0.8); //constructor for line: TLine( x1, y1, x2, y2)
    line->SetLineColor(kBlack);
    line->SetLineStyle(2);
    line->SetLineWidth(2);
    line->Draw("same");
  }
//   if (tot_max>0.) {  //<- condition never met in gastof_simple
//     TLine* line = new TLine(tot_max, 0., tot_max, h_tot_all_channels->GetMaximum()*0.8);
//     line->SetLineColor(kBlack);
//     line->SetLineStyle(2);
//     line->SetLineWidth(2);
//     line->Draw("same");
  }

//.........这里部分代码省略.........
   // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
   Float_t var1, var2;
   reader->AddVariable( "var1", &var1 );
   reader->AddVariable( "var2", &var2 );

   // Spectator variables declared in the training have to be added to the reader, too
   Float_t spec1,spec2;
   reader->AddSpectator( "spec1:=var1*2",  &spec1 );
   reader->AddSpectator( "spec2:=var1*3",  &spec2 );

   // --- Book the MVA methods

   TString dir    = "weights/";
   TString prefix = "TMVARegression";

   // Book method(s)
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = it->first + " method";
         TString weightfile = dir + prefix + "_" + TString(it->first) + ".weights.xml";
         reader->BookMVA( methodName, weightfile ); 
      }
   }
   
   // Book output histograms
   TH1* hists[100];
   Int_t nhists = -1;
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      TH1* h = new TH1F( it->first.c_str(), TString(it->first) + " method", 100, -100, 600 );
      if (it->second) hists[++nhists] = h;
   }
   nhists++;
   
   // Prepare input tree (this must be replaced by your data source)
   // in this example, there is a toy tree with signal and one with background events
   // we'll later on use only the "signal" events for the test in this example.
   //   
   TFile *input(0);
   TString fname = "./tmva_reg_example.root";
   if (!gSystem->AccessPathName( fname )) {
      input = TFile::Open( fname ); // check if file in local directory exists
   } 
   else { 
      input = TFile::Open( "http://root.cern.ch/files/tmva_reg_example.root" ); // if not: download from ROOT server
   }
   
   if (!input) {
      std::cout << "ERROR: could not open data file" << std::endl;
      exit(1);
   }
   std::cout << "--- TMVARegressionApp        : Using input file: " << input->GetName() << std::endl;

   // --- Event loop

   // Prepare the tree
   // - here the variable names have to corresponds to your tree
   // - you can use the same variables as above which is slightly faster,
   //   but of course you can use different ones and copy the values inside the event loop
   //
   TTree* theTree = (TTree*)input->Get("TreeR");
   std::cout << "--- Select signal sample" << std::endl;
   theTree->SetBranchAddress( "var1", &var1 );
   theTree->SetBranchAddress( "var2", &var2 );

   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {

      if (ievt%1000 == 0) {
         std::cout << "--- ... Processing event: " << ievt << std::endl;
      }

      theTree->GetEntry(ievt);

      // Retrieve the MVA target values (regression outputs) and fill into histograms
      // NOTE: EvaluateRegression(..) returns a vector for multi-target regression

      for (Int_t ih=0; ih<nhists; ih++) {
         TString title = hists[ih]->GetTitle();
         Float_t val = (reader->EvaluateRegression( title ))[0];
         hists[ih]->Fill( val );    
      }
   }
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();

   // --- Write histograms

   TFile *target  = new TFile( "TMVARegApp.root","RECREATE" );
   for (Int_t ih=0; ih<nhists; ih++) hists[ih]->Write();
   target->Close();

   std::cout << "--- Created root file: \"" << target->GetName() 
             << "\" containing the MVA output histograms" << std::endl;
  
   delete reader;
    
   std::cout << "==> TMVARegressionApplication is done!" << std::endl << std::endl;
}

//.........这里部分代码省略.........
   t->SetBranchAddress("b1_E", &b1_E, &b_b1_E);
   t->SetBranchAddress("b2_dR", &b2_dR, &b_b2_dR);
   t->SetBranchAddress("b2_csv", &b2_csv, &b_b2_csv);
   t->SetBranchAddress("b2_px", &b2_px, &b_b2_px);
   t->SetBranchAddress("b2_py", &b2_py, &b_b2_py);
   t->SetBranchAddress("b2_pz", &b2_pz, &b_b2_pz);
   t->SetBranchAddress("b2_E", &b2_E, &b_b2_E);
   t->SetBranchAddress("tauvis1_dR", &tauvis1_dR, &b_tauvis1_dR);
   t->SetBranchAddress("tauvis1_px", &tauvis1_px, &b_tauvis1_px);
   t->SetBranchAddress("tauvis1_py", &tauvis1_py, &b_tauvis1_py);
   t->SetBranchAddress("tauvis1_pz", &tauvis1_pz, &b_tauvis1_pz);
   t->SetBranchAddress("tauvis1_E", &tauvis1_E, &b_tauvis1_E);
   t->SetBranchAddress("tauvis2_dR", &tauvis2_dR, &b_tauvis2_dR);
   t->SetBranchAddress("tauvis2_px", &tauvis2_px, &b_tauvis2_px);
   t->SetBranchAddress("tauvis2_py", &tauvis2_py, &b_tauvis2_py);
   t->SetBranchAddress("tauvis2_pz", &tauvis2_pz, &b_tauvis2_pz);
   t->SetBranchAddress("tauvis2_E", &tauvis2_E, &b_tauvis2_E);
   t->SetBranchAddress("met_pt", &met_pt, &b_met_pt);
   t->SetBranchAddress("met_px", &met_px, &b_met_px);
   t->SetBranchAddress("met_py", &met_py, &b_met_py);
   t->SetBranchAddress("met_cov00", &met_cov00, &b_met_cov00);
   t->SetBranchAddress("met_cov10", &met_cov10, &b_met_cov10);
   t->SetBranchAddress("met_cov01", &met_cov01, &b_met_cov01);
   t->SetBranchAddress("met_cov11", &met_cov11, &b_met_cov11);

  //define the testd hypotheses
  std::vector<Int_t> hypo_mh1;
  hypo_mh1.push_back(125);
  std::vector<Int_t> hypo_mh2;
  hypo_mh2.push_back(125);
  
  // event loop
  for (int i = 0; i < max_nevent; i++) {
    t->GetEntry(i);
    
    //use only btaged jets and check for correct gen match
    if (b1_csv<0.681 || b2_csv<0.681 || njets < 2 || b1_dR > 0.2 || b2_dR > 0.2 || tauvis1_dR > 0.1 || tauvis2_dR > 0.1)   continue;
    
    //define input vectors
    TLorentzVector b1      = TLorentzVector(b1_px,b1_py,b1_pz,b1_E);
    TLorentzVector b2      = TLorentzVector(b2_px,b2_py,b2_pz,b2_E);
    TLorentzVector tau1vis = TLorentzVector(tauvis1_px,tauvis1_py,tauvis1_pz,tauvis1_E);
    TLorentzVector tau2vis = TLorentzVector(tauvis2_px,tauvis2_py,tauvis2_pz,tauvis2_E);
    
    TLorentzVector ptmiss  = TLorentzVector(met_px,met_py,0,met_pt);
    TMatrixD metcov(2,2);
    metcov(0,0)=met_cov00;
    metcov(1,0)=met_cov10;
    metcov(0,1)=met_cov01;    
    metcov(1,1)=met_cov11;
    
    //intance of fitter master class
    HHKinFitMaster kinFits = HHKinFitMaster(&b1,&b2,&tau1vis,&tau2vis);
    kinFits.setAdvancedBalance(&ptmiss,metcov);
    //kinFits.setSimpleBalance(ptmiss.Pt(),10); //alternative which uses only the absolute value of ptmiss in the fit
    kinFits.addMh1Hypothesis(hypo_mh1);
    kinFits.addMh2Hypothesis(hypo_mh2);
    kinFits.doFullFit();

    //obtain results from different hypotheses
    Double_t chi2_best = kinFits.getBestChi2FullFit();
    Double_t mh_best = kinFits.getBestMHFullFit();
    std::pair<Int_t, Int_t> bestHypo = kinFits.getBestHypoFullFit();
    std::map< std::pair<Int_t, Int_t>, Double_t> fit_results_chi2 = kinFits.getChi2FullFit();
    std::map< std::pair<Int_t, Int_t>, Double_t> fit_results_fitprob = kinFits.getFitProbFullFit();
    std::map< std::pair<Int_t, Int_t>, Double_t> fit_results_mH = kinFits.getMHFullFit();

//.........这里部分代码省略.........
		double ptZjj[99];
		double lep[99];
		int    nCands;
		int    nXjets[99];
		double region[99];

		// event weight
		double weight=-1;

		tree->SetBranchAddress("nLooseMu", &nLooseMu);
		tree->SetBranchAddress("nLooseEle",&nLooseEle);
		tree->SetBranchAddress("nsubjetbtagL",nsubjetbtagL);
		tree->SetBranchAddress("nsubjetbtagM",nsubjetbtagM);
		tree->SetBranchAddress("nsubjetbtagT",nsubjetbtagT);
		tree->SetBranchAddress("nfatjetbtagL",nfatjetbtagL); 
		tree->SetBranchAddress("nfatjetbtagM",nfatjetbtagM); 
		tree->SetBranchAddress("nfatjetbtagT",nfatjetbtagT);

		tree->SetBranchAddress("ptlep1", ptlep1);
		tree->SetBranchAddress("ptZll", ptZll);
		tree->SetBranchAddress("ptZjj", ptZjj);
		tree->SetBranchAddress("met", &met);
		tree->SetBranchAddress("nXjets", nXjets);
		tree->SetBranchAddress("lep", lep);
		tree->SetBranchAddress("nCands", &nCands);
		tree->SetBranchAddress("region", &region);

		tree->SetBranchAddress("weight", &weight);

		bool filled = 0;

		for(int j=0;j<entries;j++)
		{
			tree->GetEntry(j);

			//per event weight
			double actualWeight = weight*19538.85;
			if(i<ndata) actualWeight =1; //for data
			//cout<<actualWeight<<endl;

			filled = 0;
			//if(nCands==0)continue;
			for(int ivec=0;ivec<nCands;ivec++){

				if(filled==0)
				{
					//cout<<nXjets[ivec]<<endl;
					//make cuts
					if((nLooseEle+nLooseMu)!=1)continue;
					//if(met<80)continue;
					if(lep[ivec]!=1.)continue;//mu
					if(nXjets[ivec]!=1)continue;
					if(region[ivec]!=1)continue;
					total = total + actualWeight;

					if(nsubjetbtagL[ivec]>=1&&nfatjetbtagL[ivec]>=1)pass1=pass1+actualWeight;
					if(nsubjetbtagM[ivec]>=1&&nfatjetbtagL[ivec]>=1)pass2=pass2+actualWeight;
					if(nsubjetbtagT[ivec]>=1&&nfatjetbtagL[ivec]>=1)pass3=pass3+actualWeight;
					if(nsubjetbtagL[ivec]>=2&&nfatjetbtagL[ivec]>=1)pass4=pass4+actualWeight;
					if(nsubjetbtagM[ivec]>=2&&nfatjetbtagL[ivec]>=1)pass5=pass5+actualWeight;
					if(nsubjetbtagT[ivec]>=2&&nfatjetbtagL[ivec]>=1)pass6=pass6+actualWeight;
					if(nsubjetbtagL[ivec]>=1&&nfatjetbtagM[ivec]>=1)pass7=pass7+actualWeight;
					if(nsubjetbtagM[ivec]>=1&&nfatjetbtagM[ivec]>=1)pass8=pass8+actualWeight;
					if(nsubjetbtagT[ivec]>=1&&nfatjetbtagM[ivec]>=1)pass9=pass9+actualWeight;
					if(nsubjetbtagL[ivec]>=2&&nfatjetbtagM[ivec]>=1)pass10=pass10+actualWeight;
					if(nsubjetbtagM[ivec]>=2&&nfatjetbtagM[ivec]>=1)pass11=pass11+actualWeight;

void MergeRootfile(std::map<std::pair<std::string, std::string>, TObject*>& outputMap, std::vector<std::pair<std::string, TObject*> >& outputVec, TDirectory *target, TFile *source)
{
    using namespace std;
    string path(target->GetPath());
    path = path.substr(path.find(":") + 1);

    source->cd(path.c_str());
    TDirectory *current_sourcedir = gDirectory;
    //gain time, do not add the objects in the list in memory
    Bool_t status = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);

    // loop over all keys in this directory
    TIter nextkey( current_sourcedir->GetListOfKeys() );
    TKey *key, *oldkey = 0;
    while(key = (TKey*)nextkey())
    {
        //keep only the highest cycle number for each key
        if (oldkey && !strcmp(oldkey->GetName(), key->GetName())) continue;
        oldkey = key;

        // read object from source file
        source->cd(path.c_str());
        TObject *obj = key->ReadObj();

        if(obj->IsA()->InheritsFrom(TH1::Class()))
        {
            if(verbosity >= 4) 
            {
                printf("| Found TH1: %s\n", obj->GetName());
                fflush(stdout);
            }
            string path(target->GetPath());
            pair<string, string> okey(path.substr(path.find(':') + 2), obj->GetName());
            //cout << okey.first << "\t" << okey.second << endl;
            if(outputMap.find(okey) == outputMap.end())
            {
                outputVec.push_back(make_pair(path.substr(path.find(':') + 2), obj));
                outputMap[okey] = obj;
            }
            else
            {
                ((TH1*)outputMap[okey])->Add((TH1*)obj);
                ((TH1*)obj)->Delete();
            }
        }
        else if(obj->IsA()->InheritsFrom(TTree::Class()))
        {
            string path(target->GetPath());
            if(verbosity >= 4) 
            {
                printf("| Found Tree: %s\n", obj->GetName());
                fflush(stdout);
            }
            pair<string, string> okey(path.substr(path.find(':') + 2), obj->GetName());
            if(outputMap.find(okey) == outputMap.end())
            {
                string fname(okey.first);
                fname = fname + "_" + obj->GetName() + "tmpfile";
                for(size_t pos; (pos = fname.find('/')) != size_t(-1); ) fname[pos] = '_';
                tmpFiles.push_back(make_pair(fname, new TFile(fname.c_str(), "RECREATE")));
                TTree* tree = ((TTree*)obj)->CloneTree();
                ((TTree*)obj)->GetListOfClones()->Remove(tree);
                ((TTree*)obj)->ResetBranchAddresses();
                tree->ResetBranchAddresses();
                outputVec.push_back(make_pair(path.substr(path.find(':') + 2), (TObject*)tree));
                outputMap[okey] = (TObject*)tree;
            }
            else 
            {
                TTree* tm = (TTree*)outputMap[okey];
                TTree* ts = (TTree*)obj;
                tm->CopyAddresses(ts);
                for(int i = 0; i < ts->GetEntries(); i++)
                {
                    ts->GetEntry(i);
                    tm->Fill();
                }
                ts->ResetBranchAddresses();
                if (tm->GetTreeIndex()) tm->GetTreeIndex()->Append(ts->GetTreeIndex(), kTRUE); 
                ((TTree*)obj)->Delete();
            }
        }
        else if(obj->IsA()->InheritsFrom(TDirectory::Class()))
        {
            if(verbosity >= 3) 
            {
                printf("Hadding Directory: %s\n", ((TDirectory*)obj)->GetPath());
                fflush(stdout);
            }
            string path(((TDirectory*)obj)->GetPath());
            pair<string, string> okey(path.substr(path.find(':') + 2), " -------- ");
            if(outputMap.find(okey) == outputMap.end())
            {
                outputVec.push_back(make_pair(path.substr(path.find(':') + 2), (TDirectory*)0));
                outputMap[okey] = 0;
            }
            MergeRootfile(outputMap, outputVec, (TDirectory*)obj, source);
        }
        else
//.........这里部分代码省略.........

double TrimEventContent(Int_t iRapBin = 1,
		      Int_t iPTBin = 1,
		      Double_t fracL = 0.5, Double_t nSigma = 2., 
		      Int_t nUpsState=0,//[0]... 1S, [1]... 2S, [2]... 3S
		      bool UpsMC=false,
		      bool f_BG_zero=false,
		      bool ProjectLSBdata=false,
		      bool ProjectRSBdata=false,
		      bool CombineSignalPeaks=false,
		      bool Y1Sto2S_SB=false,
		      bool LeftSided=false,
		      bool RightSided=false,
		      bool MassScan=false,
		      bool adjustOverlapBorders=true
		      ){

  printf("\n\n\nfracL = %1.3f, nSigma = %1.1f, iState = %d, rap %d, pT %d\n", fracL, nSigma, nUpsState, iRapBin, iPTBin);

  Char_t name[100], title[100];
  Char_t fileNameIn[100];
  sprintf(fileNameIn, "tmpFiles/data_Ups_rap%d_pT%d.root", iRapBin, iPTBin);
  //==============================
  //read inputs from input file:
  TFile *fIn = new TFile(fileNameIn);
  TLorentzVector *lepP;
  TLorentzVector *lepN;
  TTree *treeIn = (TTree *) gDirectory->Get("selectedData");
  if(gDirectory->Get("selectedData")==NULL){
    printf("\n\n\nskip processing this bin.\n\n\n");
    return -999.;
  }

  TH2D *hBG_cosThetaPhiLR[onia::kNbFrames][2];
  for(int iFrame = 0; iFrame < onia::kNbFrames; iFrame++){
    sprintf(name, "hBG_cosThetaPhi_%s_L", onia::frameLabel[iFrame]);
    hBG_cosThetaPhiLR[iFrame][L] = (TH2D *) gDirectory->Get(name);
    sprintf(name, "hBG_cosThetaPhi_%s_R", onia::frameLabel[iFrame]);
    hBG_cosThetaPhiLR[iFrame][R] = (TH2D *) gDirectory->Get(name);
  }
  //==============================

  //definition of output variables 
  Char_t fileNameOut[100];
  sprintf(fileNameOut, "AllStates_%1.2fSigma_FracLSB%dPercent/data_%dSUps_rap%d_pT%d.root", nSigma, int(fracL*100), nUpsState+1, iRapBin, iPTBin);
  TFile *fOut = new TFile(fileNameOut, "RECREATE");
  gStyle->SetPadRightMargin(0.2);
  TTree *treeOut = treeIn->CloneTree(0);
  // treeOut->SetName("data");
  TH2D *hBG_cosThetaPhi[onia::kNbFrames];
  // TH2D *hBG_cosThetaPhiSignal[onia::kNbFrames];
  for(int iFrame = 0; iFrame < onia::kNbFrames; iFrame++){
    // //book the histo for the signal
    // sprintf(name, "total_%s", onia::frameLabel[iFrame]);
    // sprintf(title, ";cos#theta_{%s};#phi_{%s} [deg]", onia::frameLabel[iFrame], onia::frameLabel[iFrame]);
    // hBG_cosThetaPhiSignal[iFrame] = new TH2D(name, title, onia::kNbBinsCosT, onia::cosTMin, onia::cosTMax, 
    // 					  onia::kNbBinsPhiPol, onia::phiPolMin, onia::phiPolMax);
    // hBG_cosThetaPhiSignal[iFrame]->Sumw2();
    //copy the L and R sideband histos into one output BG histogram
    hBG_cosThetaPhiLR[iFrame][L]->Scale(fracL/hBG_cosThetaPhiLR[iFrame][L]->Integral());
    hBG_cosThetaPhiLR[iFrame][R]->Scale((1.-fracL)/hBG_cosThetaPhiLR[iFrame][R]->Integral());
    sprintf(name, "background_costhphi%s", onia::frameLabel[iFrame]);
    hBG_cosThetaPhi[iFrame] = (TH2D *) hBG_cosThetaPhiLR[iFrame][L]->Clone(name);
    hBG_cosThetaPhi[iFrame]->Add(hBG_cosThetaPhiLR[iFrame][R]);
  }

  //==========================================================
  //reading fit parameters to establish signal mass window
  //as well as the L and R sideband window for the 3D BG histo
  //==========================================================
  fIn->cd();
  TTree *treeFitPar = (TTree *) gDirectory->Get("massFitParameters");
  TF1 *fUps[kNbSpecies], *fBG = 0;
  fUps[0] = 0, fUps[1] = 0, fUps[2] = 0;
  treeFitPar->SetBranchAddress("fUps1S", &fUps[0]);
  treeFitPar->SetBranchAddress("fUps2S", &fUps[1]);
  treeFitPar->SetBranchAddress("fUps3S", &fUps[2]);
  treeFitPar->SetBranchAddress("fBG", &fBG);
  treeFitPar->LoadTree(0);
  treeFitPar->GetEntry(0);


  Double_t mass[kNbSpecies], sigma[kNbSpecies];
  for(int iState = 0; iState < kNbSpecies; iState++){
    mass[iState] = fUps[iState]->GetParameter(1);
    sigma[iState] = fUps[iState]->GetParameter(2);
  }
  printf("1S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS1S], sigma[UPS1S]);
  printf("2S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS2S], sigma[UPS2S]);
  printf("3S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS3S], sigma[UPS3S]);
  Double_t poleMass = mass[nUpsState], massMin, massMax;
  massMin = poleMass - nSigma*sigma[nUpsState];
  massMax = poleMass + nSigma*sigma[nUpsState];

  if(LeftSided){
	  massMin = poleMass - nSigma*sigma[nUpsState];
	  massMax = poleMass;
  }

  if(RightSided){
	  massMin = poleMass;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	UpsilonTree_allsign->Branch("muPlusPt", &muPlusPt, "muPlusPt/F");
	UpsilonTree_allsign->Branch("muMinusPt", &muMinusPt, "muMinusPt/F");
	UpsilonTree_allsign->Branch("muPlusEta", &muPlusEta, "muPlusEta/F");
	UpsilonTree_allsign->Branch("muMinusEta", &muMinusEta, "muMinusEta/F");
	UpsilonTree_allsign->Branch("muPlusPhi", &muPlusPhi, "muPlusPhi/F");
	UpsilonTree_allsign->Branch("muMinusPhi", &muMinusPhi, "muMinusPhi/F");
	//	UpsilonTree_allsign->Branch("theta_HX", &theta_HX, "theta_HX/F");
	//	UpsilonTree_allsign->Branch("phi_HX", &phi_HX, "phi_HX/F");
	//	UpsilonTree_allsign->Branch("theta_CS", &theta_CS, "theta_CS/F");
	//	UpsilonTree_allsign->Branch("phi_CS", &phi_CS, "phi_CS/F");

	UpsilonTree_trkRot->Branch("Centrality",    &Centrality,    "Centrality/I");
	UpsilonTree_trkRot->Branch("HLTriggers",    &HLTriggers,    "HLTriggers/I");
	UpsilonTree_trkRot->Branch("QQtrig",   &QQtrig,    "QQtrig/I");
	UpsilonTree_trkRot->Branch("QQsign",   &QQsign,    "QQsign/I");
	UpsilonTree_trkRot->Branch("weight",   &weight,    "weight/F");
	UpsilonTree_trkRot->Branch("weight2",   &weight2,    "weight2/F");
	UpsilonTree_trkRot->Branch("vProb",    &vProb,     "vProb/F");
	UpsilonTree_trkRot->Branch("eventNb",    &eventNb,     "eventNb/I");
	UpsilonTree_trkRot->Branch("runNb",    &runNb,     "runNb/I");
	UpsilonTree_trkRot->Branch("invariantMass", &invariantMass, "invariantMass/F");
	UpsilonTree_trkRot->Branch("upsPt", &upsPt, "upsPt/F");
	UpsilonTree_trkRot->Branch("upsEta", &upsEta, "upsEta/F");
	UpsilonTree_trkRot->Branch("upsRapidity", &upsRapidity, "upsRapidity/F");
	UpsilonTree_trkRot->Branch("muPlusPt", &muPlusPt, "muPlusPt/F");
	UpsilonTree_trkRot->Branch("muMinusPt", &muMinusPt, "muMinusPt/F");
	UpsilonTree_trkRot->Branch("muPlusEta", &muPlusEta, "muPlusEta/F");
	UpsilonTree_trkRot->Branch("muMinusEta", &muMinusEta, "muMinusEta/F");
	UpsilonTree_trkRot->Branch("muPlusPhi", &muPlusPhi, "muPlusPhi/F");
	UpsilonTree_trkRot->Branch("muMinusPhi", &muMinusPhi, "muMinusPhi/F");


	for (int i=0; i<nentries; i++) {
		t->GetEntry(i);
		//nReco_QQ=0;
		if (i%1000==0) cout<<i<<endl;
		nPlusMu=0;
		nMinusMu=0;
		for(int iMu = 0; iMu < Reco_mu_size; iMu++){
			if (Reco_mu_charge[iMu] == 1) nPlusMu++;
			else nMinusMu++;
			TLorentzVector *Reco_mu = (TLorentzVector *) Reco_mu_4mom->At(iMu);
			muPt[iMu]=Reco_mu->Pt();
			muEta[iMu]=Reco_mu->Eta();
			muPhi[iMu]=Reco_mu->Phi();
		}
		MuTree->Fill();

		for(int iQQ = 0; iQQ < Reco_QQ_size; iQQ++){
			//eventNb=eventNb[iQQ];
			//runNb=runNb[iQQ];
			vProb = Reco_QQ_VtxProb[iQQ];
			QQtrig = Reco_QQ_trig[iQQ];
			QQsign = Reco_QQ_sign[iQQ];
			if (Reco_QQ_sign[iQQ] == 0) { 
				weight = 1;
				weight2 = 1;
			}
			else {
				weight = -1;
				float likesign_comb = (float)nPlusMu*(nPlusMu-1.0)/2.0+(float)nMinusMu*(nMinusMu-1.0)/2.0;
				float unlikesign_bkgd_comb = (float)nPlusMu*nMinusMu - (nPlusMu>nMinusMu?nMinusMu:nPlusMu);
				weight2 = -1.0 * unlikesign_bkgd_comb/likesign_comb;
			}

			TLorentzVector *Reco_QQ = (TLorentzVector *) Reco_QQ_4mom->At(iQQ);

void kanaelec::Loop(int wda, int wsp, const char *outfilename)
{
  if (fChain == 0) return;
  Long64_t nentries = fChain->GetEntries();
  // Out Put File Here
  TFile fresults = TFile(outfilename,"RECREATE");
  TTree *newtree = fChain->CloneTree();
  
  Float_t fit_mu_px=0,   fit_mu_py =0,   fit_mu_pz=0,   fit_mu_e=0;
  Float_t fit_nv_px=0,   fit_nv_py =0,   fit_nv_pz=0,   fit_nv_e=0;
  Float_t fit_aj_px=0,   fit_aj_py =0,   fit_aj_pz=0,   fit_aj_e=0;
  Float_t fit_bj_px=0,   fit_bj_py =0,   fit_bj_pz=0,   fit_bj_e=0;
  Float_t fit_mlvjj=0,   fit_chi2  =999;  
  Int_t   fit_NDF  =999, fit_status=999, gdevtt   =0;

  TBranch *branch_mu_px = newtree->Branch("fit_el_px", &fit_mu_px,  "fit_el_px/F");
  TBranch *branch_mu_py = newtree->Branch("fit_el_py", &fit_mu_py,  "fit_el_py/F");
  TBranch *branch_mu_pz = newtree->Branch("fit_el_pz", &fit_mu_pz,  "fit_el_pz/F");
  TBranch *branch_mu_e  = newtree->Branch("fit_el_e",  &fit_mu_e,   "fit_el_e/F");
  	   
  TBranch *branch_nv_px = newtree->Branch("fit_nv_px", &fit_nv_px,  "fit_nv_px/F");
  TBranch *branch_nv_py = newtree->Branch("fit_nv_py", &fit_nv_py,  "fit_nv_py/F");
  TBranch *branch_nv_pz = newtree->Branch("fit_nv_pz", &fit_nv_pz,  "fit_nv_pz/F");
  TBranch *branch_nv_e  = newtree->Branch("fit_nv_e",  &fit_nv_e,   "fit_nv_e/F");
  	   
  TBranch *branch_aj_px = newtree->Branch("fit_aj_px", &fit_aj_px,  "fit_aj_px/F");
  TBranch *branch_aj_py = newtree->Branch("fit_aj_py", &fit_aj_py,  "fit_aj_py/F");
  TBranch *branch_aj_pz = newtree->Branch("fit_aj_pz", &fit_aj_pz,  "fit_aj_pz/F");
  TBranch *branch_aj_e  = newtree->Branch("fit_aj_e",  &fit_aj_e,   "fit_aj_e/F");
  	   
  TBranch *branch_bj_px = newtree->Branch("fit_bj_px", &fit_bj_px,  "fit_bj_px/F");
  TBranch *branch_bj_py = newtree->Branch("fit_bj_py", &fit_bj_py,  "fit_bj_py/F");
  TBranch *branch_bj_pz = newtree->Branch("fit_bj_pz", &fit_bj_pz,  "fit_bj_pz/F");
  TBranch *branch_bj_e  = newtree->Branch("fit_bj_e",  &fit_bj_e,   "fit_bj_e/F");
  
  TBranch *branch_mlvjj = newtree->Branch("fit_mlvjj", &fit_mlvjj,  "fit_mlvjj/F");
  TBranch *branch_chi2  = newtree->Branch("fit_chi2",  &fit_chi2,   "fit_chi2/F");
  TBranch *branch_NDF   = newtree->Branch("fit_NDF",   &fit_NDF,    "fit_NDF/I");
  TBranch *branch_status= newtree->Branch("fit_status",&fit_status, "fit_status/I");
  TBranch *branch_gdevtt= newtree->Branch("gdevtt",    &gdevtt,     "gdevtt/I");

  Float_t fi2_mu_px=0,   fi2_mu_py =0,   fi2_mu_pz=0,   fi2_mu_e=0;
  Float_t fi2_nv_px=0,   fi2_nv_py =0,   fi2_nv_pz=0,   fi2_nv_e=0;
  Float_t fi2_aj_px=0,   fi2_aj_py =0,   fi2_aj_pz=0,   fi2_aj_e=0;
  Float_t fi2_bj_px=0,   fi2_bj_py =0,   fi2_bj_pz=0,   fi2_bj_e=0;
  Float_t fi2_mlvjj=0,   fi2_chi2  =999;  
  Int_t   fi2_NDF  =999, fi2_status=999;

  TBranch *branc2_mu_px = newtree->Branch("fi2_el_px", &fi2_mu_px,  "fi2_el_px/F");
  TBranch *branc2_mu_py = newtree->Branch("fi2_el_py", &fi2_mu_py,  "fi2_el_py/F");
  TBranch *branc2_mu_pz = newtree->Branch("fi2_el_pz", &fi2_mu_pz,  "fi2_el_pz/F");
  TBranch *branc2_mu_e  = newtree->Branch("fi2_el_e",  &fi2_mu_e,   "fi2_el_e/F");
  	   
  TBranch *branc2_nv_px = newtree->Branch("fi2_nv_px", &fi2_nv_px,  "fi2_nv_px/F");
  TBranch *branc2_nv_py = newtree->Branch("fi2_nv_py", &fi2_nv_py,  "fi2_nv_py/F");
  TBranch *branc2_nv_pz = newtree->Branch("fi2_nv_pz", &fi2_nv_pz,  "fi2_nv_pz/F");
  TBranch *branc2_nv_e  = newtree->Branch("fi2_nv_e",  &fi2_nv_e,   "fi2_nv_e/F");
  	   
  TBranch *branc2_aj_px = newtree->Branch("fi2_aj_px", &fi2_aj_px,  "fi2_aj_px/F");
  TBranch *branc2_aj_py = newtree->Branch("fi2_aj_py", &fi2_aj_py,  "fi2_aj_py/F");
  TBranch *branc2_aj_pz = newtree->Branch("fi2_aj_pz", &fi2_aj_pz,  "fi2_aj_pz/F");
  TBranch *branc2_aj_e  = newtree->Branch("fi2_aj_e",  &fi2_aj_e,   "fi2_aj_e/F");
  	   
  TBranch *branc2_bj_px = newtree->Branch("fi2_bj_px", &fi2_bj_px,  "fi2_bj_px/F");
  TBranch *branc2_bj_py = newtree->Branch("fi2_bj_py", &fi2_bj_py,  "fi2_bj_py/F");
  TBranch *branc2_bj_pz = newtree->Branch("fi2_bj_pz", &fi2_bj_pz,  "fi2_bj_pz/F");
  TBranch *branc2_bj_e  = newtree->Branch("fi2_bj_e",  &fi2_bj_e,   "fi2_bj_e/F");
  
  TBranch *branc2_mlvjj = newtree->Branch("fi2_mlvjj", &fi2_mlvjj,  "fi2_mlvjj/F");
  TBranch *branc2_chi2  = newtree->Branch("fi2_chi2",  &fi2_chi2,   "fi2_chi2/F");
  TBranch *branc2_NDF   = newtree->Branch("fi2_NDF",   &fi2_NDF,    "fi2_NDF/I");
  TBranch *branc2_status= newtree->Branch("fi2_status",&fi2_status, "fi2_status/I");

  // Loop over all events
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    //Long64_t ientry = LoadTree(jentry);
    //if (ientry < 0) break;
    nb = newtree->GetEntry(jentry);   nbytes += nb;
    // Cut variable definitions
    double jess    = 1.00; // control the jet energy scale
    //double muoniso = (W_muon_trackiso+W_muon_hcaliso+W_muon_ecaliso-event_RhoForLeptonIsolation*3.141592653589*0.09)/W_muon_pt;
    double dijetpt = sqrt(JetPFCor_Pt[0]*JetPFCor_Pt[0]+
			  JetPFCor_Pt[1]*JetPFCor_Pt[1]+
			  2*JetPFCor_Pt[0]*JetPFCor_Pt[1]*cos(JetPFCor_Phi[0]-JetPFCor_Phi[1]));

    // Save variable initialization
    fit_mu_px = 0; fit_mu_py = 0; fit_mu_pz = 0;  fit_mu_e  = 0; 
    fit_nv_px = 0; fit_nv_py = 0; fit_nv_pz = 0;  fit_nv_e  = 0; 
    fit_aj_px = 0; fit_aj_py = 0; fit_aj_pz = 0;  fit_aj_e  = 0; 
    fit_bj_px = 0; fit_bj_py = 0; fit_bj_pz = 0;  fit_bj_e  = 0; 
    fit_mlvjj = 0; fit_chi2  =999;fit_NDF   =999; fit_status=999;

    gdevtt    = 0;

    fi2_mu_px = 0; fi2_mu_py = 0; fi2_mu_pz = 0;  fi2_mu_e  = 0; 
    fi2_nv_px = 0; fi2_nv_py = 0; fi2_nv_pz = 0;  fi2_nv_e  = 0; 
    fi2_aj_px = 0; fi2_aj_py = 0; fi2_aj_pz = 0;  fi2_aj_e  = 0; 
    fi2_bj_px = 0; fi2_bj_py = 0; fi2_bj_pz = 0;  fi2_bj_e  = 0; 
    fi2_mlvjj = 0; fi2_chi2  =999;fi2_NDF   =999; fi2_status=999;
//.........这里部分代码省略.........