C++ TEventList类代码示例

////////////////////////
// Clear the event list
SEXP clear(SEXP eventList)
{
  TEventList* el = checkForEventListWrapper(eventList);
  
  el->Clear();
  
  return R_NilValue;
}

///////////////////////
// Get the number of entries in the event list
SEXP nEntriesEventList(SEXP eventList)
{
  TEventList* el = checkForEventListWrapper(eventList);
  
  SEXP n = NEW_INTEGER(1);
  INTEGER(n)[0] = el->GetN();
  
  return n;
}

///////////////////////
// Get the title of the event list
SEXP getTitle(SEXP eventList)
{
  TEventList* el = checkForEventListWrapper(eventList);
  
  SEXP name = NEW_CHARACTER(1);
  SET_STRING_ELT( name, 0, mkChar(el->GetTitle()) );
  
  return name;
}

/////////////////////////
// Put entries into the event list
SEXP putIntoEventList(SEXP eventList, SEXP entryNums)
{
  
  TEventList* el = checkForEventListWrapper(eventList);
  
  for ( unsigned int i = 0; i < GET_LENGTH(entryNums); ++i ) {
    el->Enter( INTEGER(entryNums)[i] );
  }

  return R_NilValue;
}

////////////////
// Get the event list
SEXP RootChainManager::getEventListName()
{
  TEventList* ev = m_chain->GetEventList();
	
  if (! ev ) return R_NilValue;
	
  SEXP r = NEW_CHARACTER(1);
  SET_STRING_ELT(r, 0, mkChar(ev->GetName()) );
	
  return r;
}

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 a logical vector if entries are in the event list
SEXP isInEventList(SEXP eventList, SEXP entryNums)
{
  
  TEventList* el = checkForEventListWrapper(eventList);
  
  SEXP l;
  PROTECT(l = NEW_LOGICAL( GET_LENGTH(entryNums) ) );
  
  for ( unsigned int i = 0; i < GET_LENGTH(entryNums); ++i ) {
    LOGICAL(l)[i] = el->Contains( INTEGER(entryNums)[i] ) == 1;
  }
  
  UNPROTECT(1);
  
  return l;
}

void treeTest(  const TString fList = "temp.list",
                const TString fOutName = "fOut.root",
                Long64_t nEvents = 9999
                ) 
{

	gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
    gROOT->LoadMacro("makeEventList.C");
	loadSharedLibraries();

    TFile* fOut = new TFile(fOutName.Data(),"RECREATE");

    TH1D* hZdc = new TH1D("hZdc","hZdc", 50, 0.5, 5000.5);

	//----------------- Read data using StMuDstMaker ---------------//
	StChain* muChain = new StChain("StChain");
	StMuDstMaker* muMaker = new StMuDstMaker(0,0,"",fList.Data(),"st:MuDst.root", 100, "MuDst");

    muChain->Init();
    Float_t zdcCut = 1623;

    TEventList* eventList = makeEventList(muMaker->chain(), zdcCut);
    muMaker->SetEventList(eventList);
    Long64_t nList = eventList->GetN();
    if(nList < nEvents) {nEvents = nList;}

    for(Int_t i = 0; i <= (nEvents-1); i++)
    {
        muChain->Make();
        Float_t zdcWest = muMaker->muDst()->event()->zdcTriggerDetector().adc(0);
        Float_t zdcEast = muMaker->muDst()->event()->zdcTriggerDetector().adc(4);

        hZdc->Fill(zdcWest);
        hZdc->Fill(zdcEast);


    }

    muChain->Finish();
    fOut->Write();
    fOut->Close();

    delete muChain;
    delete eventList;

}

//////////////////////////
// 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;
}

    int CEntityParent::lTriggerEvent(lua_State* L)
    {
        //first param is userdata
        int nEventId = luaL_checkint(L, 2);
        CHECK_UNSIGNED(nEventId, 2);

        //参数个数
        int nParam = lua_gettop(L) - 2;

        if (nParam >= 0)
        {
            CEventDispatcher * p = GetWorld()->GetEventDispatcher();
            TEventList* l = p->TriggerEvent(GetId(), nEventId);
            if (l)
            {
                TEventList::iterator iter1 = l->begin();
                map<TENTITYID, string> tmap;
                for (;iter1 != l->end(); iter1++)
                {
                    tmap.insert(make_pair(iter1->first, iter1->second));
                }

                map<TENTITYID, string>::iterator iter2 = tmap.begin();
                for (;iter2 != tmap.end(); ++iter2)
                {
                    CEntityParent* pe = GetWorld()->GetEntity(iter2->first);
                    if(pe)
                    {
                        int n = EntityMethodCall(L, pe, iter2->second.c_str(), nParam, 0);
                        lua_pop(L, n);
                        //LogDebug("lTriggerEvent", "%d", l->size());
                    }
                }
            }
        }

        return 0;
    }

//.........这里部分代码省略.........
		<< "trigger.HLT_PFHT350_PFMET100_v3==1 || trigger.HLT_PFHT350_PFMET100_v4==1 || trigger.HLT_PFHT350_PFMET100_v5==1 || "
		<< "trigger.HLT_PFHT350_PFMET100_v6==1 || trigger.HLT_PFHT350_PFMET100_v7==1 || trigger.HLT_PFNoPUHT350_PFMET100_v1==1 || "
		<< "trigger.HLT_PFNoPUHT350_PFMET100_v3==1 || trigger.HLT_PFNoPUHT350_PFMET100_v4==1 ) ) )&&TOBTECTagger<=8&&ExtraBeamHaloFilter==0)";

	TString trigger = triggerStream.str().c_str();
	load(samples.Data());

	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;

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();
}

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

void FitSignals(TTree * treeB, TCut cut, Int_t max){
  AliSignalProcesor proc;
  TF1 * f1 = proc.GetAsymGauss();
  TTreeSRedirector cstream("FitSignal.root");
  TFile *f = cstream.GetFile();

  char lname[100];
  sprintf(lname,"Fit%s", cut.GetTitle());
  TEventList *list = new TEventList(lname,lname);
  sprintf(lname,">>Fit%s", cut.GetTitle());
  treeB->Draw(lname,cut);
  treeB->SetEventList(list);
  Int_t nFits=0;
  for (Int_t ievent=0; ievent<list->GetN(); ievent++){
    if (nFits>max) break;
    if (nFits%50==0) printf("%d\n",nFits);
    char ename[100];
    sprintf(ename,"Fit%d", ievent);
    Double_t nsample = treeB->Draw("Graph.fY-Mean09:Graph.fX","","",1,ievent);
    Double_t * signal  = treeB->GetV1();
    Double_t * time  = treeB->GetV2();
    Double_t maxpos =0;
    Double_t max = 0;
    for (Int_t ipos = 0; ipos<nsample; ipos++){
      if (signal[ipos]>max){
	max    = signal[ipos];
	maxpos = ipos;
      }
    }

    Int_t first = TMath::Max(maxpos-10,0.);
    Int_t last  = TMath::Min(maxpos+60, nsample);
    //
    f->cd();
    TH1F his(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his.SetBinContent(ipos+1,signal[ipos+first]);
    }
    treeB->Draw("Sector:Row:Pad","","",1,ievent);
    Double_t sector = treeB->GetV1()[0];
    Double_t row    = treeB->GetV2()[0];
    Double_t pad    = treeB->GetV3()[0];
    //    TGraph  graph(last-first,&time[first],&signal[first]);
    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
    //    TH1F * his = (TH1F*)graph.GetHistogram();
    his.Fit(f1,"q");
    his.Write(ename);
    gPad->Clear();
    his.Draw();
    gPad->Update();
    Double_t params[6];
    for (Int_t ipar=0; ipar<6; ipar++) params[ipar] = f1->GetParameters()[ipar];
    Double_t chi2 = TFitter::GetFitter()->Chisquare(6,params);
    TMatrixD cov(6,6);
    cov.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());
    //
    // tail cancellation
    //
    Double_t x0[1000];
    Double_t x1[1000];
    Double_t x2[1000];
    for (Int_t ipos=0; ipos<last-first; ipos++){
      x0[ipos] = signal[ipos+first];
    }
    proc.TailCancelationALTRO1(x0,x1,0.85*0.339,0.09,last-first);
    proc.TailCancelationALTRO1(x1,x2,0.85,0.789,last-first);
    //
    sprintf(ename,"Cancel1_%d", ievent);
    TH1F his1(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his1.SetBinContent(ipos+1,x1[ipos]);
    }
    his1.Write(ename);
    sprintf(ename,"Cancel2_%d", ievent);
    TH1F his2(ename,ename,last-first,first,last);
    for (Int_t ipos=0; ipos<last-first; ipos++){
      his2.SetBinContent(ipos+1,x1[ipos]);
    }
    f1->SetParameters(0.75*max,maxpos,1.1,0.8,0.25,0.2);
    his2.Fit(f1,"q");
    his2.Write(ename);
    Double_t params2[6];
    for (Int_t ipar=0; ipar<6; ipar++) params2[ipar] = f1->GetParameters()[ipar];
    Double_t chi22 = TFitter::GetFitter()->Chisquare(6,params2);    
    TMatrixD cov2(6,6);
    cov2.SetMatrixArray(TFitter::GetFitter()->GetCovarianceMatrix());

    TGraph gr0(last-first, &time[first],x0);
    TGraph gr1(last-first, &time[first],x1);
    TGraph gr2(last-first, &time[first],x2);
    //
    cstream<<"Fit"<<
      "Sector="<<sector<<
      "Row="<<row<<
      "Pad="<<pad<<
      "First="<<first<<
      "Max="<<max<<
      "MaxPos="<<maxpos<<
      "chi2="<<chi2<<
      "chi22="<<chi22<<
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  T->SetAlias("HCALIN_MIP", "Energy_Sum_Hadron_HCALIN<0.5");
  T->SetAlias("Energy_Sum_Hadron_HCALOUT",
              "12./9.50430e+00*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");

  T->SetAlias("MIP_Count_Col2",
              "Sum$( abs( TOWER_RAW_CEMC.get_energy() )>20 && abs( TOWER_RAW_CEMC.get_energy() )<200 && TOWER_CALIB_CEMC.get_column() == 2 )");
  T->SetAlias("Pedestal_Count_AllCEMC",
              "Sum$( abs( TOWER_RAW_CEMC.get_energy() )<20)");

  //
  TCut event_sel = "1*1";

  if (plot_all)
  {
    //      event_sel = "1*1";
    //      cuts = "_all_event";
    //      event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL";
    //      cuts = "_Valid_HODO";
    event_sel =
        "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO";
    cuts = "_Valid_HODO_Trigger_VETO";
  }
  else
  {
    if (0)
    {  // energy selection
      event_sel = Form("(beam_MTNRG_GeV == %f)", beam_momentum_selection);
      cuts = Form("_%.0fGeV", beam_momentum_selection);

      cout << "Build event selection of " << (const char *) event_sel
           << endl;

      T->Draw(">>EventListRunCut", event_sel);
      TEventList *elist = gDirectory->GetObjectChecked("EventListRunCut",
                                                       "TEventList");
      cout << elist->GetN() << " / " << T->GetEntriesFast()
           << " events selected" << endl;
      T->SetEventList(elist);
    }

    //      event_sel = "1*1";
    //      cuts = "_all_event";
    //      event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL";
    //      cuts = "_Valid_HODO";
    event_sel =
        "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO";
    cuts = "_Valid_HODO_Trigger_VETO";

    //      event_sel =
    //          event_sel
    //              + "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO  && (C2_Sum_e<50)";
    //      cuts = cuts + "_Valid_HODO_Trigger_VETO_C2_Sum_Hadron";

    //      event_sel =
    //          event_sel
    //              + "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && CEMC_MIP && (C2_Sum_e<100)";
    //      cuts = cuts + "_Valid_HODO_Trigger_VETO_CEMC_MIP_C2_Sum_Hadron";

    //      event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && (MIP_Count_Col2 == 8) && (Pedestal_Count_AllCEMC >= 64 - 8)";
    //      cuts = "_Valid_HODO_MIP_Col2_PedestalOther";

    //      event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && Average_HODO_VERTICAL>1.5 && Average_HODO_VERTICAL<4.5 && Average_HODO_HORIZONTAL>3.5 && Average_HODO_HORIZONTAL<6.5";
    //      cuts = "_Valid_HODO_center_col1_row2";

    //      event_sel =
    //          "C2_Inner_e>100 && Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && Average_HODO_HORIZONTAL>1.5 && Average_HODO_HORIZONTAL<4.5 && Average_HODO_VERTICAL>3.5 && Average_HODO_VERTICAL<6.5";