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

void NCIdeogram::addTree(TTree *tree, const char *mu, const char *sigma)
{
	TLeaf *leafMu = tree->FindLeaf(mu);
	TLeaf *leafSigma = tree->FindLeaf(sigma);
	TEventList *elist = tree->GetEventList();
	Long64_t j,ix,nbr;
	
	nbr = elist ? elist->GetN() : tree->GetEntries();
	for (j = 0; j < nbr; j++) {
		ix = elist ? elist->GetEntry(j) : j;
		tree->GetEntry(ix);
		addPoint(leafMu->GetValue(), leafSigma->GetValue());
	}
} // addTree()

//////////////////////////
// Return all of the events as a big vector
SEXP getEntries(SEXP eventList)
{
  TEventList* el = checkForEventListWrapper(eventList);
  
  unsigned int l = el->GetN();
  
  SEXP e;
  PROTECT( e = NEW_INTEGER(l) );
  
  for ( unsigned int i = 0; i < l; ++i) {
    INTEGER(e)[i] = el->GetEntry(i);
  }
  
  UNPROTECT(1);
  
  return e;
}

void listMatchedEvents(std::string fullPath, bool matchesonly = false)
{
    // Essentially a sort of diff for reco and gen tree
    // May get slow for large trees as the search for the matching entry in the second tree
    // traverses the whole tree. This is because for merged trees the order is not guaranteed
    const int fVerbose(1);
    // Open file
    TFile *f = TFile::Open(fullPath.c_str());
    if (f==0)
    {
	cout << "File " << fullPath << " not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0)
	cout << "Succesfully opened file " << fullPath << endl;
    // Get TTree with GenEvents
    TTree* tgen = (TTree*) f->Get("genevents");
    if(fVerbose>0) cout << "Got TTree with " << tgen->GetEntries() << " entries" << endl;
    // Do a cut, if needed
    //tgen->Draw(">>lst","ptmu1>3&&ptmu2>3&&TMath::Abs(etamu1)<2.5&&TMath::Abs(etamu2)<2.5");
    tgen->Draw(">>lst","");
    TEventList *lst;
    lst = (TEventList*)gDirectory->Get("lst");
    tgen->SetEventList(lst);
    if(fVerbose>0) cout << " After cuts: " << lst->GetN() << " entries" << endl;

    // Get TTree with iRecoEvents
    TTree* treco = (TTree*) f->Get("events");
    if(fVerbose>0) cout << "Got TTree with " << treco->GetEntries() << " entries" << endl;

    // set branch addresses
    int genrun, genls, genevt;
    tgen->SetBranchAddress("run",&genrun);
    tgen->SetBranchAddress("LS",&genls);
    tgen->SetBranchAddress("event",&genevt);

    int recorun, recols, recoevt;
    treco->SetBranchAddress("run",&recorun);
    treco->SetBranchAddress("LS",&recols);
    treco->SetBranchAddress("event",&recoevt);

    int mcmatch;
    treco->SetBranchAddress("isMCmatch",&mcmatch);

    double ptmu1, ptmu2, ptpr, ptpi;
    tgen->SetBranchAddress("ptmu1",&ptmu1);
    tgen->SetBranchAddress("ptmu2",&ptmu2);
    tgen->SetBranchAddress("ptpr",&ptpr);
    tgen->SetBranchAddress("ptpi",&ptpi);
    double etamu1, etamu2, etapr, etapi;
    tgen->SetBranchAddress("etamu1",&etamu1);
    tgen->SetBranchAddress("etamu2",&etamu2);
    tgen->SetBranchAddress("etapr",&etapr);
    tgen->SetBranchAddress("etapi",&etapi);
    double phimu1, phimu2, phipr, phipi;
    tgen->SetBranchAddress("phimu1",&phimu1);
    tgen->SetBranchAddress("phimu2",&phimu2);
    tgen->SetBranchAddress("phipr",&phipr);
    tgen->SetBranchAddress("phipi",&phipi);

    cout << "#      run         LS      event     genidx    recoidx" << endl;
    cout << "---------- ---------- ---------- ---------- ----------" << endl;

    TLorentzVector tlvmu1, tlvmu2, tlvpr, tlvpi;
    tmph1 = new TH1F ("tmph1","tmphisto",20,3.096,3.098);
    tmph2 = new TH1F ("tmph2","tmphisto",20,1.112,1.120);
    tmph3 = new TH1F ("tmph3","tmphisto",20,5.61,5.64);
    tmph4 = new TH2F ("tmph4","tmphisto",20,0,5,20,-2,8);

    // outer loop on gen events
    for (int i=0; i!=lst->GetN(); i++)
    {
	tgen->GetEntry(lst->GetEntry(i));
	tlvmu1.SetPtEtaPhiM(ptmu1,etamu1,phimu1,0.105658367);
	tlvmu2.SetPtEtaPhiM(ptmu2,etamu2,phimu2,0.105658367);
	tlvpr.SetPtEtaPhiM( ptpr, etapr, phipr ,0.938272013);
	tlvpi.SetPtEtaPhiM( ptpi, etapi, phipi ,0.13957018);
	tmph1->Fill((tlvmu1+tlvmu2).M());
	tmph2->Fill((tlvpr+tlvpi).M());
	tmph3->Fill((tlvmu1+tlvmu2+tlvpr+tlvpi).M());
	tmph4->Fill((tlvpr+tlvpi).P(),tlvpr.P()-tlvpi.P());
	//cout << (tlvmu1+tlvmu2).M() << " " << (tlvpr+tlvpi).M() << " " << (tlvmu1+tlvmu2+tlvpr+tlvpi).M() << endl;
	if (!matchesonly)
	    cout << setw(10) << genrun << " " << setw(10) << genls << " " << setw(10) << genevt << " " << setw(10) << lst->GetEntry(i) << " ";
	// inner loop on reco events
	int nfound(0);
	for (int j=0; j!=treco->GetEntries(); j++)
	{
	    treco->GetEntry(j);
	    if(genrun==recorun && genls==recols && genevt==recoevt)
	    {
		nfound++;
		if (matchesonly && nfound==1)
		    cout << setw(10) << genrun << " " << setw(10) << genls << " " << setw(10) << genevt << " " << setw(10) << lst->GetEntry(i) << " ";
		if(nfound>1) cout << "                                 " << setw(10) << lst->GetEntry(i) << " ";
		cout << setw(10) << j << " match: " << mcmatch << endl;
	    }
	}
	if (!matchesonly && nfound==0) cout << "         -" << endl;
    }
//.........这里部分代码省略.........

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

	int counterAll = 0;

  	for(size_t i = 0; i < fSamples.size(); ++i){

	    string sampletype = (string)fSamples[i].type;
	    if(sampletype!="data") continue;//run only over data - extendable to MC 

	    Double_t sample_weight = fSamples[i].xsection * fSamples[i].kfact * fSamples[i].lumi / (fSamples[i].nevents*fSamples[i].PU_avg_weight);
	    if(fVerbose>2) cout << "ZnunuNumbers: looping over " << fSamples[i].name << " added in " << sampletype << endl;
	    if(fVerbose>2) cout << "              sample has weight " << sample_weight << " and " << fSamples[i].tree->GetEntries() << " entries" << endl; 
	    if(fVerbose>2 && fSamples[i].tree->GetEntries()==0) cout << "skip sample, has no entries" << endl;
            if(fSamples[i].tree->GetEntries()==0) continue;

	    MT2tree* fMT2tree = new MT2tree();
	    fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
	    Long64_t nentries =  fSamples[i].tree->GetEntries();
	    Long64_t nbytes = 0, nb = 0;
	    int nev =0;
        
	    TString myCuts = cuts;// + "&&" + basecuts;
        
	    if( fSamples[i].type=="data") myCuts += " && " + trigger; //cuts to be aplied only on data
        
   	    cout << "Cuts for Flow: " << myCuts << endl;
   	    fSamples[i].tree->Draw(">>selList", myCuts);
        
	    TEventList *myEvtList = (TEventList*)gDirectory->Get("selList");
	    fSamples[i].tree->SetEventList(myEvtList);
	    int counter=0;
	    cout << "Filtering done, size=" <<myEvtList->GetN()  << endl;
	    if(myEvtList->GetSize()==0) continue;
	*fLogStream << "Sample " << fSamples[i].name << endl;
	//run over all selected events
        while(myEvtList->GetEntry(counter++) !=-1){	
		int jentry = myEvtList->GetEntry(counter-1);
            
		nb =  fSamples[i].tree->GetEntry(jentry);   nbytes += nb;
		fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
		
		if ( fVerbose>2 && counter % 5000 == 0  )  cout << "+++ Proccessing event " << counter << endl;

		string sHT;//HT region
		if(fMT2tree->misc.HT<450.)      sHT = "_HTge0";
		else if(fMT2tree->misc.HT<750.) sHT = "_lowHT";
		else if(fMT2tree->misc.HT<1200.)sHT = "_mediumHT";
		else                            sHT = "_highHT";

		string ssignal;//topological region
		if(fMT2tree->NJetsIDLoose40 == 2 &&                                  fMT2tree->NBJets40CSVM == 0) ssignal = "_2j0b";
		if(fMT2tree->NJetsIDLoose40 == 2 &&                                  fMT2tree->NBJets40CSVM >= 1) ssignal = "_2j1to2b";
		if(fMT2tree->NJetsIDLoose40 >= 3 && fMT2tree->NJetsIDLoose40 <= 5 && fMT2tree->NBJets40CSVM == 0) ssignal = "_3to5j0b";
		if(fMT2tree->NJetsIDLoose40 >= 3 && fMT2tree->NJetsIDLoose40 <= 5 && fMT2tree->NBJets40CSVM == 1) ssignal = "_3to5j1b";
		if(fMT2tree->NJetsIDLoose40 >= 3 && fMT2tree->NJetsIDLoose40 <= 5 && fMT2tree->NBJets40CSVM == 2) ssignal = "_3to5j2b";
		if(                                                                  fMT2tree->NBJets40CSVM >= 3) ssignal = "_3b";
		if(fMT2tree->NJetsIDLoose40 >= 6 &&                                  fMT2tree->NBJets40CSVM == 0) ssignal = "_6j0b";
		if(fMT2tree->NJetsIDLoose40 >= 6 &&                                  fMT2tree->NBJets40CSVM == 1) ssignal = "_6j1b";
		if(fMT2tree->NJetsIDLoose40 >= 6 &&                                  fMT2tree->NBJets40CSVM == 2) ssignal = "_6j2b";

		string hh = sHT + ssignal;
		histos[(string)"MT2"        + hh]->Fill(fMT2tree->misc.MT2);//fill histogram
		if(fMT2tree->misc.MT2>histos[(string)"MT2"        + hh]->GetBinLowEdge(1)){
			if(sHT=="_lowHT" && (fMT2tree->trigger.HLT_PFHT350_PFMET100_v3==1||fMT2tree->trigger.HLT_PFHT350_PFMET100_v4==1||fMT2tree->trigger.HLT_PFHT350_PFMET100_v5==1||fMT2tree->trigger.HLT_PFHT350_PFMET100_v6==1||fMT2tree->trigger.HLT_PFHT350_PFMET100_v7==1 || fMT2tree->trigger.HLT_PFNoPUHT350_PFMET100_v1==1 || fMT2tree->trigger.HLT_PFNoPUHT350_PFMET100_v3==1 || fMT2tree->trigger.HLT_PFNoPUHT350_PFMET100_v4==1 ) ) *fLogStream << "HTMHT ";
			else if(sHT=="_lowHT") *fLogStream << "MET   ";
			else *fLogStream << "JetHT ";
			*fLogStream << int(fMT2tree->misc.Run) << ":" << int(fMT2tree->misc.LumiSection) << ":" << int(fMT2tree->misc.Event) << endl;//file log file (including trigger stream
			++counterAll;
		}
	}//while
	delete fMT2tree;
	delete fSamples[i].tree;

	}//for samples

	cout << "add overflow to last bin" << endl;
	for(map<string,TH1D*>::iterator h=histos.begin(); h!=histos.end();++h){
		h->second->SetBinContent(h->second->GetNbinsX(),
					 h->second->GetBinContent(h->second->GetNbinsX()  )+ 
					 h->second->GetBinContent(h->second->GetNbinsX()+1)  );
		h->second->SetBinError(  h->second->GetNbinsX(),
					 sqrt(h->second->GetBinError(h->second->GetNbinsX()  )*
					      h->second->GetBinError(h->second->GetNbinsX()  )+
					      h->second->GetBinError(h->second->GetNbinsX()+1)*
					      h->second->GetBinError(h->second->GetNbinsX()+1)  ));
	}

	cout << "Saving." << endl;
    	TFile *fsavefile = new TFile(fOutDir + outputname,"RECREATE");
	fsavefile->cd();
	for(map<string,TH1D*>::iterator h=histos.begin(); h!=histos.end();++h){
		h->second->Write();
	}
	fsavefile->Close();
	cout << "Saved histograms in " << fOutDir << outputname << endl;

	TString logname ="eventlistMT2inclusive.log"; 
	cout << counterAll << " events in " << logname.Data() << endl;
	ofstream f_log (logname.Data(), ios::trunc);
	f_log << fLogStream->str();
}

void Plot(Char_t *filen = "jsf.root")
{
  //gROOT->Reset();
  gSystem->Load("libS4Utils.so");
  gSystem->Load("libAnlib.so");

  Int_t nZoneX = 4;
  Int_t nZoneY = 3;

  TCanvas *c1 = new TCanvas("c1","",0,0, 300*nZoneX, 300*nZoneY);
  c1->SetHighLightColor(5);
  c1->SetFillColor(19);
  c1->Divide(nZoneX,nZoneY);

  gStyle->SetOptFit();

  cerr << filen << endl;
  TFile   *filep = new TFile(filen);
  TNtupleD *tup   = (TNtupleD *)gROOT->FindObject("hEvt");

  //--
  // Preselection
  //--
#if 0
  TString sel("evis>0.");
#else
  TString sel("pt<20.&&abs(pl)<30.&&ycut>0.0015");
#endif
  tup->Draw(">>elist",sel.Data(),"goff");

  TEventList *elist = static_cast<TEventList*>(gROOT->FindObject("elist"));
  Int_t nlist = elist->GetN();
  cerr << "Nevent = " << nlist << endl;

  //--
  // Loop over preselected events
  //--
  Double_t kToDeg = 180./TMath::Pi();
  for (Int_t i=0; i<nlist; i++) {
#if 0
        Int_t event = elist->GetEntry(i);
        cerr << "i = " << i << " event = " << event << endl;
#else
        Int_t event = i;
#endif
        Double_t ntracks;
        Double_t evis;
        Double_t pt;
        Double_t pl;
        Double_t ycut;
        Double_t chi2;
        Double_t mh;
        Double_t mw1;
        Double_t mw2;
        Double_t mt1;
        Double_t mt2;
        Double_t mtt;

        tup->SetBranchAddress("ntracks",&ntracks);
        tup->SetBranchAddress("evis",&evis);
        tup->SetBranchAddress("pt",&pt);
        tup->SetBranchAddress("pl",&pl);
        tup->SetBranchAddress("ycut",&ycut);
        tup->SetBranchAddress("chi2",&chi2);
        tup->SetBranchAddress("mh",&mh);
        tup->SetBranchAddress("mw1",&mw1);
        tup->SetBranchAddress("mw2",&mw2);
        tup->SetBranchAddress("mt1",&mt1);
        tup->SetBranchAddress("mt2",&mt2);
        tup->SetBranchAddress("mtt",&mtt);

        tup->GetEntry(event);

        hMhMw->Fill(mh,  mw1, 1.);
        hMwMw->Fill(mw1, mw2, 1.);
        hMtMt->Fill(mt1, mt2, 1.);
	hEvis->Fill(evis, 1.);
	hPtPl->Fill(pt, pl, 1.);
        hMttMh->Fill(mtt,  mh, 1.);
  }
  Int_t id = 0;
  id++; c1->cd(id); hStat->Draw();
  id++; c1->cd(id); hMhMw->Draw();
  id++; c1->cd(id); hMhMw->ProjectionX()->Draw();
  id++; c1->cd(id); hMhMw->ProjectionY()->Draw();
  id++; c1->cd(id); hEvis->Draw();
  id++; c1->cd(id); hMwMw->Draw();
  id++; c1->cd(id); hMwMw->ProjectionX()->Draw();
  id++; c1->cd(id); hMwMw->ProjectionY()->Draw();
  id++; c1->cd(id); hMtMt->Draw();
  id++; c1->cd(id); hMtMt->ProjectionX()->Draw();
  id++; c1->cd(id); hMtMt->ProjectionY()->Draw();
  id++; c1->cd(id); hMttMh->ProjectionX()->Draw();
}

//call this macro via root -l -b -q AddBranchEventsMT2tree.C++
//This code takes a text file containing run:lumisection:event:taggervalue (tagger value is a new filter you want to add to your tree),
//an old MT2tree file that does not contain that tagger
// and creates a new MT2tree file, which is a copy of the old one and contains additionally that tagger.
void AddBranchEventsMT2tree() {
// Example of Root macro to copy a subset of a Tree to a new Tree
// Only selected entries are copied to the new Tree.
// The input file has been generated by the program in $ROOTSYS/test/Event
// with   Event 1000 1 99 1
//Author: Rene Brun
   
//modified: Hannsjoerg Weber, 28/08/2013

   //make cleaning
   // get eventlist to clean
   // vector contains <run, < <lumisection, event>, taggervalue > >
   vector<pair<pair<int,pair<int,int> >,float> > rls; rls.clear();
   char buffer[200];
   ifstream filterdat("/shome/haweber/AODFiles_MT2SR/taggerlist/ControlRegionSinglePhotonPart1.dat");
   while( filterdat.getline(buffer, 200, '\n') ){
	int rrun(-1), lls(-1), eevent(-1), d1(-1); float ttagger(-1);
	sscanf(buffer, "*\t%d\t*\t%d\t*\t%d\t*\t%d\t*\t%f", &d1, &eevent, &lls, &rrun, &ttagger);
	pair<int,int> t1(lls,eevent);
	pair<int,pair<int,int> > t2(rrun,t1);
        pair<pair<int,pair<int,int> >, float> t3(t2,ttagger);
	rls.push_back(t3);
   }
cout << "Events for filtering " << rls.size() << endl;

   // get the old rootfile including events that should be tagged and define new event file including tagger variable
   TString oldfilename = "/shome/haweber/MT2Analysis/MT2trees/MT2_V02-03-02/20130914_8TeV_1g_removed/lowHT/SinglePhoton-Run2012C-PromptReco-v2-2.root";
   TString newfilename = "/shome/haweber/MT2Analysis/MT2trees/MT2_V02-03-02/20130914_8TeV_1g_removed/lowHT/SinglePhoton-Run2012C-PromptReco-v2-2_Filter.root";

	// if you want to apply an additional event selection, define cuts here
	std::ostringstream cutStream;
	std::ostringstream cutStreamBase;
	cutStream << " " 
		<< "NTausIDLoose3Hits+NMuons+NEles==0"                   << "&&"
		<< "misc.Jet0Pass==1&&misc.Jet1Pass==1"                  << "&&"
		<< "misc.Vectorsumpt<70&& misc.MinMetJetDPhi4Pt40 >0.3";
		cutStream << "&&((misc.MET>200&&misc.HT<=750&&misc.HT>=450&&misc.MT2>200)||(misc.HT>750&&misc.MET>30&&misc.MT2>=100))";
	cutStreamBase << " " 
      << "misc.PassJet40ID ==1"                             << "&&"
      << "(misc.HBHENoiseFlag == 0 || misc.ProcessID==10)"  << "&&" // not there for fastsim (i.e. signal)
      << "misc.CSCTightHaloIDFlag == 0"                     << "&&"
      << "misc.trackingFailureFlag==0"                      << "&&"
      << "misc.eeBadScFlag==0"                              << "&&"
      << "misc.EcalDeadCellTriggerPrimitiveFlag==0"         << "&&"
      << "misc.TrackingManyStripClusFlag==0"                << "&&"
      << "misc.TrackingTooManyStripClusFlag==0"             << "&&"
      << "misc.TrackingLogErrorTooManyClustersFlag==0"      << "&&"
      << "misc.CrazyHCAL==0";

	TString cuts = cutStream.str().c_str();
	TString basecuts = cutStreamBase.str().c_str();
//	TString myCuts = cuts + "&&" + basecuts;// if cuts should be applied DO NOT comment this line
	TString myCuts = "";// if no cuts should be applied DO NOT comment this line

   //Get old file, old tree and set top branch address
   TFile *oldfile = new TFile(oldfilename.Data());
   TTree *oldtree = (TTree*)oldfile->Get("MassTree");

   MT2tree* fMT2tree = new MT2tree();
   oldtree->SetBranchAddress("MT2tree", &fMT2tree);
   Long64_t nentries =  oldtree->GetEntries();
   Long64_t nbytes = 0, nb = 0;
   int nev =0;

   //Create a new file + a clone of old tree in new file
   Float_t TOBTECTagger(-1);
   Float_t HOTagger(-1);
   Bool_t ExtraBeamHaloFilter(false);
   TFile *newfile = new TFile(newfilename.Data(),"RECREATE");
   TTree *newtree = oldtree->CloneTree(0);//clone all branches, but no event
   TBranch *newBranch  = newtree->Branch("TOBTECTagger", &TOBTECTagger, "TOBTECTagger/F");
   TBranch *newBranch2 = newtree->Branch("ExtraBeamHaloFilter", &ExtraBeamHaloFilter, "ExtraBeamHaloFilter/O");
   TBranch *newBranch3 = newtree->Branch("HOTagger", &HOTagger, "HOTagger/F");

   oldtree->Draw(">>selList", myCuts);
   TEventList *myEvtList = (TEventList*)gDirectory->Get("selList");
   oldtree->SetEventList(myEvtList);
   int counter=0;
   int numEvt = myEvtList->GetN();
   int printEvt = 100000;
   if(myEvtList->GetN()<1000000) printEvt = myEvtList->GetN()/25;
   cout << "Filtering done, size=" <<myEvtList->GetN()  << endl;
   while(myEvtList->GetEntry(counter++) !=-1){
	int jentry = myEvtList->GetEntry(counter-1);
	nb =  oldtree->GetEntry(jentry);   nbytes += nb;
	oldtree->SetBranchAddress("MT2tree", &fMT2tree);
	if(counter%printEvt==0) cout << "Process event " << counter << endl;	

	//make cleaning
	bool keep = true;
	int p=-1;
	for(int nn=0; nn<rls.size();++nn){
		if((rls[nn].first).first!=fMT2tree->misc.Run) continue; // --> run over matching 
		if(((rls[nn].first).second).first!=fMT2tree->misc.LumiSection) continue; // --> LS over matching 
		if(fMT2tree->misc.Event>0){
			if(((rls[nn].first).second).second==fMT2tree->misc.Event) { 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	    Double_t sample_weight = fSamples[i].xsection * fSamples[i].kfact * fSamples[i].lumi / (fSamples[i].nevents);
	    if(fVerbose>2) cout << "B-Efficiency: looping over " << fSamples[i].name << endl;
	    if(fVerbose>2) cout << "              sample has weight " << sample_weight << " and " << fSamples[i].tree->GetEntries() << " entries" << endl; 
	    if(fVerbose>2 && fSamples[i].tree->GetEntries()==0) cout << "skip sample, has no entries" << endl;
            if(fSamples[i].tree->GetEntries()==0) continue;
	    MT2tree* fMT2tree = new MT2tree();
	    fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
	    Long64_t nentries =  fSamples[i].tree->GetEntries();
	    Long64_t nbytes = 0, nb = 0;
	    int nev =0;

	    TString myCuts = cuts;
	    if(HTskim)  myCuts = myCuts + " && misc.HT>750";
	    if(METskim) myCuts = myCuts + " && misc.HT<=750 && misc.MET>200";//METskim excluding HTskim

	/*
	    if( fSamples[i].type=="data" && fSamples[i].sname=="EE-Data") { myCuts += " && " + triggerEE; }//cuts to be aplied only on data
	    else if( fSamples[i].type=="data" && fSamples[i].sname=="EMu-Data") { myCuts += " && " + triggerEMu; }//cuts to be aplied only on data
	    else if( fSamples[i].type=="data" && fSamples[i].sname=="MuMu-Data") { myCuts += " && " + triggerMuMu; } //cuts to be aplied only on data
	    else if(fSamples[i].type=="data") {cout << "data not usuable" << " type " << fSamples[i].type << " Sname " << fSamples[i].sname << endl; continue; }//not 
	*/
	    if( fSamples[i].type=="data") { myCuts += " && " + trigger; }

   	    cout << "Cuts for Flow: " << myCuts << endl;
   	    fSamples[i].tree->Draw(">>selList", myCuts);

	    TEventList *myEvtList = (TEventList*)gDirectory->Get("selList");

	    fSamples[i].tree->SetEventList(myEvtList);

	    int counter=0;
	    cout << "Filtering done, size=" <<myEvtList->GetN()  << endl;
	    if(myEvtList->GetSize()==0) continue;
        while(myEvtList->GetEntry(counter++) !=-1){	
      		int jentry = myEvtList->GetEntry(counter-1);
            
            nb =  fSamples[i].tree->GetEntry(jentry);   nbytes += nb;
            fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
            
            if ( fVerbose>2 && counter % 50000 == 0 )  cout << "+++ Proccessing event " << counter << endl;


	    //if((fSamples[i].sname!="VVV"&&fSamples[i].sname!="TTbarV" && fSamples[i].type!="susy"&&sampletype!="Stop")  && fMT2tree->misc.HBHENoiseFlag!=0     ) continue;//test
	    //if((fSamples[i].sname!="VVV"&&fSamples[i].sname!="TTbarV" && fSamples[i].type!="susy"&&sampletype!="Stop")  && fMT2tree->misc.hcalLaserEventFlag!=0) continue;//test
            Double_t weight = sample_weight;
            if (!fMT2tree->misc.isData ) weight = weight * fMT2tree->pileUp.Weight; // pile-up reweighting for MC 
 /*           Bool_t recoedee   = false;// exact 2 ele, 0 muo
            Bool_t recoedemu  = false;// exact 1 ele, 1 muo
            Bool_t recoedmumu = false;// exact 0 ele, 2 muo
		//change this - might have some inefficiencies - three leptons with one not passing IDMedium==1
		//20-20 selection --> third lepton is 10 GeV 
            if(fMT2tree->NEles>=2 &&fMT2tree->ele[0].lv.Pt()>20&&fMT2tree->ele[1].lv.Pt()>20&&fMT2tree->ele[2].lv.Pt()<20&&(fMT2tree->NMuons==0||fMT2tree->muo[0].lv.Pt()<10) && fMT2tree->ele[0].IDMedium==1 && fMT2tree->ele[1].IDMedium==1) recoedee   = true;
            if(fMT2tree->NMuons>=2&&fMT2tree->muo[0].lv.Pt()>20&&fMT2tree->muo[1].lv.Pt()>20&&fMT2tree->muo[2].lv.Pt()<10&&(fMT2tree->NEles ==0||fMT2tree->ele[0].lv.Pt()<10)) recoedmumu = true;
            if(fMT2tree->NMuons>=1&&fMT2tree->muo[0].lv.Pt()>20&&fMT2tree->muo[1].lv.Pt()<10&&fMT2tree->NEles>=1 &&fMT2tree->ele[0].lv.Pt()>20&&fMT2tree->ele[1].lv.Pt()<10 && fMT2tree->ele[0].IDMedium==1) recoedemu  = true;//XXX change emu to 20/20 selection, left ee,mumu at 20/10 for now
	    //id cuts?

            Bool_t recoedosee   = false;// opposite sign
            Bool_t recoedosemu  = false;// opposite sign
            Bool_t recoedosmumu = false;// opposite sign
	    if(recoedee   && (fMT2tree->ele[0].Charge)*(fMT2tree->ele[1].Charge)==(-1)) recoedosee   = true;
	    if(recoedemu  && (fMT2tree->ele[0].Charge)*(fMT2tree->muo[0].Charge)==(-1)) recoedosemu  = true;
	    if(recoedmumu && (fMT2tree->muo[0].Charge)*(fMT2tree->muo[1].Charge)==(-1)) recoedosmumu = true;
            Bool_t recoedeenZ   = false;// off Z
            Bool_t recoedemunZ  = false;// off Z -> is not requirement, so this variable should be useless
            Bool_t recoedmumunZ = false;// off Z
	    if(recoedee   && ((fMT2tree->ele[0].lv + fMT2tree->ele[1].lv).M()<81 || (fMT2tree->ele[0].lv + fMT2tree->ele[1].lv).M()>101) ) recoedeenZ   = true;

//this code create the so-called 'PrintOuts' as done with the MT2tree::PrintOut() function
void PrintOutEvents(){

	bool savetofile = true;//save the printout to the file PrintOut.log, default = true
	bool onlyData = true;  //run only over data, default = true

//	TString  samples = "/shome/haweber/MT2Analysis_8TeV/Code/MT2AnalysisCode/RootMacros/samples/samples_MET_filter.dat";//only dummy
//	TString  samples = "/shome/haweber/MT2Analysis_8TeV/Code/MT2AnalysisCode/RootMacros/samples/samples_HT_filter.dat";//only dummy
	TString samples = "samples/samples_HTandMET_filter.dat";

	//event selection for which you want to have the printout - complicated stuff can also be done later
	std::ostringstream cutStream;
	std::ostringstream cutStreamBase;
	cutStream << " " 
	//	<< "NBJets40CSVM >= 3"		<< "&&"
		<< "NJetsIDLoose40 >= 2"		<< "&&"
		<< "NEles+NMuons+NTausIDLoose3Hits  ==0"<< "&&"
	//	<< "misc.MT2>=250"                      << "&&"
		<< "misc.Jet0Pass ==1"                      << "&&"
		<< "misc.Jet1Pass ==1"                      << "&&"
		<< "misc.Vectorsumpt < 70";
		cutStream  << "&& misc.MinMetJetDPhi4Pt40 >0.3";
//	cutStream << "&&misc.HT<750&&misc.HT>=450&&misc.MET>200";
//	cutStream << "&&misc.HT>=750&&misc.MET>30";
	cutStream << "&&((misc.HT<750&&misc.HT>=450&&misc.MET>200)||(misc.HT>=750&&misc.MET>30))";

	//Higgs
//	cutStream <<"&&NJetsIDLoose40 >= 4"<<"&&"
  //  		  <<"NBJetsCSVM>=2"<<"&&"
	//	  <<"GetSelBBMinv()>250";
	
	cutStreamBase << " " 
      << "misc.PassJet40ID ==1"                      << "&&"
      << "(misc.HBHENoiseFlag == 0 || misc.ProcessID==10)"  << "&&" // for rare SM samples
      << "misc.CSCTightHaloIDFlag == 0"             << "&&"
      << "misc.trackingFailureFlag==0"              << "&&"
      << "misc.eeBadScFlag==0"                      << "&&"
      << "misc.EcalDeadCellTriggerPrimitiveFlag==0" << "&&"
      << "misc.TrackingManyStripClusFlag==0"             << "&&"
      << "misc.TrackingTooManyStripClusFlag==0"          << "&&"
      << "misc.TrackingLogErrorTooManyClustersFlag==0"   << "&&"
      << "misc.CrazyHCAL==0";
	cutStreamBase << "&&misc.MET/misc.CaloMETRaw<=2.";
	cutStreamBase << "&&TOBTECTagger<=8&&ExtraBeamHaloFilter==0";// ONLY ON DATA
	TString cuts = cutStream.str().c_str();
	TString basecuts = cutStreamBase.str().c_str();

	load(samples.Data());

	//if have extremely many events can only select a random subset of those
	TRandom3 *random = new TRandom3(0);//somehow truly random with this seed
	vector<int> randvec; randvec.clear();
	for(int i =0; i<40;++i)
		randvec.push_back(random->Integer(7750));
	delete random;
	sort(randvec.begin(), randvec.end());
	//for(int i=0; i<randvec.size();++i) cout << "randvec["<<i<<"]="<<randvec[i] << endl;
	int counter1 = 0; int counter2 = 0;

	//don't run events twice - this vector makes that sure (in case event is in both MET and HT trigger stream)
	   vector<pair<int,pair<int,int> > > rls; rls.clear();


   	for(size_t i = 0; i < fSamples.size(); ++i){
        
   	    if(onlyData && fSamples[i].type!="data") continue;

	    MT2tree* fMT2tree = new MT2tree();
	    fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
	    Long64_t nentries =  fSamples[i].tree->GetEntries();
	    Long64_t nbytes = 0, nb = 0;
	    int nev =0;
        
	    TString myCuts = cuts + "&&" + basecuts;
        
	    //if( fSamples[i].type=="data") myCuts += " && " + trigger; //cuts to be aplied only on data
        
	    cout << "sample " << fSamples[i].name << endl;
   	    cout << "Cuts for Flow: " << myCuts << endl;
   	    fSamples[i].tree->Draw(">>selList", myCuts);
        
	    TEventList *myEvtList = (TEventList*)gDirectory->Get("selList");
	    fSamples[i].tree->SetEventList(myEvtList);
	    int counter=0;
	    cout << "Filtering done, size=" <<myEvtList->GetN()  << endl;
	    if(myEvtList->GetSize()==0) continue;
		//run over selected events
        while(myEvtList->GetEntry(counter++) !=-1){	
		int jentry = myEvtList->GetEntry(counter-1);
            
		nb =  fSamples[i].tree->GetEntry(jentry);   nbytes += nb;
		fSamples[i].tree->SetBranchAddress("MT2tree", &fMT2tree);
		//implementation of running over random subset
//		if(counter1==randvec[counter2]){
//			++counter2;
//			fMT2tree->PrintOut(savetofile);
//		}

		pair <int,int> dummy1(fMT2tree->misc.LumiSection,fMT2tree->misc.Event);
		pair <int,pair<int,int> > dummy(fMT2tree->misc.Run, dummy1);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      input = TFile::Open( fname ); // check if file in local directory exists
   
   if (!input) {
      std::cout << "ERROR: could not open data file" << std::endl;
      exit(1);
   }
   std::cout << "--- TMVAClassificationApp    : Using input file: " << input->GetName() << std::endl;
   
   // --- Event loop

   // Prepare the event 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
   //
   std::cout << "--- Select signal sample" << std::endl;
   TTree* theTree = (TTree*)input->Get("MVA");
   //Float_t userVar1, userVar2;
   theTree->SetBranchAddress( "met",              &met );
   theTree->SetBranchAddress( "HT",               &HT );
   theTree->SetBranchAddress( "minMLB",           &minMLB );
   theTree->SetBranchAddress( "leptonJetsMETSum", &leptonJetsMETSum );


   // spectators
   theTree->SetBranchAddress( "leptonSumMass", &leptonSumMass );
   theTree->SetBranchAddress( "nJets", &nJets );
   theTree->SetBranchAddress( "nBTag", &nBTag );
   theTree->SetBranchAddress( "nLeptons", &nLeptons );
   theTree->SetBranchAddress( "isMuon1", &isMuon1 );
   theTree->SetBranchAddress( "isMuon2", &isMuon2 );
   theTree->SetBranchAddress( "isMuon3", &isMuon3 );
   theTree->SetBranchAddress( "isMuon4", &isMuon4 );


   std::vector<Float_t> vecVar(4); // vector for EvaluateMVA tests


   // CMS STATS:
   //       
   //       A little trick: loop over only selected events

   // make event list
   theTree->Draw(">>EvtList", mycut);
   TEventList * pEvtList = (TEventList *)gROOT->FindObject("EvtList");
   long int nEvents = pEvtList->GetN();

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

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

      //theTree->GetEntry(ievt);
      theTree->GetEntry(pEvtList->GetEntry(ievt));


      // CMS STATS:
      // *** HERE we assign integer branches' values to float spectator variables
      //     Otherwise our category cuts would fail
      //
      nBTagFloat=(Float_t)nBTag;
      nJetsFloat=(Float_t)nJets;
      nLeptonsFloat=(Float_t)nLeptons;
      isMuon1Float  = (Float_t)isMuon1;
      isMuon2Float  = (Float_t)isMuon2;
      isMuon3Float  = (Float_t)isMuon3;
      isMuon4Float  = (Float_t)isMuon4;



      // --- Return the MVA outputs and fill into histograms

      if (Use["BDT"          ])   histBdt    ->Fill( reader->EvaluateMVA( "BDT method"           ) );
      if (Use["Category"     ])   histCat    ->Fill( reader->EvaluateMVA( "Category method"      ) );

   }

   // Get elapsed time
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();


   // --- Write histograms

   TFile *target  = new TFile( "TMVApp.root","RECREATE" );
   if (Use["BDT"          ])   histBdt    ->Write();
   if (Use["Category"     ])   histCat    ->Write();

   target->Close();

   std::cout << "--- Created root file: \"TMVApp.root\" containing the MVA output histograms" << std::endl;
  
   delete reader;
    
   std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
}