#include <algorithm>
#include <iostream>
#include <map>
#include <vector>
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TEventList.h"
#include "TDirectory.h"
#include "TLine.h"
#include "TArrow.h"
#include "TStyle.h"
#include "TObjArray.h"
#include "TRegexp.h"

Typedefs
typedef std::vector< TString >	StrVec

typedef std::map< double, TString >	DblStrMap

typedef std::map< double, TString >::iterator	DblStrMapIt

typedef std::map< TString, double >	StrDblMap

typedef std::map< TString, int >	StrIntMap

typedef std::map< TString, int > ::iterator	StrIntMapIt

Functions
void	FCPrintInfo ()

int	FCSplitString (TString s, TString delim, StrVec &toks)

void	FCPrepareTree (TString bnames, TString ctlvar)

bool	FCFindLimits (TString var, double &min, double &max)

double	FCSfc (double S, double B)

double	FCQaVar (TString var, double qaopt=0)

void	FCInit (int numvars)

void	FCDrawVariable (TString var, int numpad, int norm=1, TString cut="")

int	findcuts (TTree *theTree=0, TString ctlvar="", TString sigcut="", TString bnames="", TString precut="", int numvars=9, double qaopt=0., double ws=1., double wb=1., int norm=1, int bins=500)

Variables
TCanvas *	c1

TEventList *	els

TEventList *	elb

TEventList *	elpre

TTree *	t

int	N0s

int	N0b

int	Ns

int	Nb

double	Wsig

double	Wbkg

TString	comtemplate

TString	comcurrent

bool	normalize

int	BINS

StrVec	bnam

DblStrMap	qamap

StrDblMap	sigeff

StrDblMap	bkgeff

Typedef Documentation

typedef std::map<double, TString> DblStrMap

Definition at line 43 of file findcuts.C.

typedef std::map<double, TString>::iterator DblStrMapIt

Definition at line 44 of file findcuts.C.

typedef std::map<TString, double> StrDblMap

Definition at line 46 of file findcuts.C.

typedef std::map<TString, int> StrIntMap

Definition at line 48 of file findcuts.C.

typedef std::map<TString, int>::iterator StrIntMapIt

Definition at line 49 of file findcuts.C.

typedef std::vector<TString> StrVec

Definition at line 41 of file findcuts.C.

Function Documentation

void FCDrawVariable	(	TString	var,
		int	numpad,
		int	norm = `1`,
		TString	cut = `""`
	)

Definition at line 390 of file findcuts.C.

References BINS, c1, comcurrent, comtemplate, cut, elb, els, FCFindLimits(), h1, h2, max(), normalize, r1, t, TString, Wbkg, Wsig, xmax, and xmin.

Referenced by findcuts().

 {
         double xmin, xmax;
         
         c1->cd(numpad);
         
         // if pad #1 (control variable), background color is set differently 
         if (numpad==1)
         {
                 gPad->SetFrameFillColor(18);
                 TString com = TString::Format(comcurrent.Data(),!normalize); // toggle normalize
                 gPad->AddExec("ex1","if (gPad->GetEvent()==kButton1Double) { cout<<endl<<\" -> Toggle normalize setting.\"<<endl;"+com+";}");
         }
         
         // for all other pads a TExec is added which when clicked applies
         // the corresponding cut 
         if (numpad!=1)
         {
                 // create the command string from the template string
                 TString com = comtemplate;
                 // isolate the part  "var>" or "var<" from the full cut string; 
                 // this is to replace a previous cut on identical variable and direction (<,>)  
                 TString cutpart = cut(0,max(cut.Index(">"),cut.Index("<"))+1);
                 TRegexp r1(cutpart+"[-]*[0-9]*.[0-9]*");
                 TString oldcut = com(r1);
                 
                 // if there was no such identical variable cut, just replace the 
                 // keyword 'NEWCUT' in template with the new cut string
                 if (oldcut=="") com.ReplaceAll("NEWCUT",cut);
                 else  // otherwise replace old cut with new cut and remove 'NEWCUT' or '&&NEWCUT' from template
                 {
                         com.ReplaceAll(oldcut,cut);
                         com.ReplaceAll("&&NEWCUT","");
                         com.ReplaceAll("NEWCUT","");
                 }
                 // finally add TExec to the current pad
                 gPad->AddExec("ex1","if (gPad->GetEvent()==kButton1Double) { cout<<endl<<\" -> Applying '"+cut+"'\"<<endl;"+com+";}");
         }
         
         // set title for histogram
         TString title = cut;
         if (title=="") title=var;
         
         // prepare line and arrow for drawing cut 
         TLine ln;
         ln.SetLineColor(6);     ln.SetLineStyle(7); ln.SetLineWidth(2);
         
         TArrow arr;
         arr.SetAngle(40); arr.SetLineWidth(2); arr.SetLineColor(6); arr.SetFillColor(6);
         
         bool drawcut = (cut!="");
         int cutdir = cut.Contains("<") ? -1 : 1 ;
         
         // determind value of cut for displaying the line and arrow
         double cutval = drawcut ? TString(cut(max(cut.Index(">"),cut.Index("<"))+1,1000)).Atof() : -999.;
                 
         if (FCFindLimits(var, xmin, xmax))
         {
                 // histograms for signal and background
                 TH1F h1("h1",title,BINS,xmin,xmax);
                 TH1F h2("h2",title,BINS,xmin,xmax);
                 // and the sum
                 TH1F hs("hs",title,BINS,xmin,xmax);
                 
                 h1.SetLineColor(4);
                 h2.SetLineColor(2);
                 hs.SetLineColor(1);
                 
                 // fill histograms
                 t->SetEventList(els);
                 t->Draw(var+">>h1");
                 t->SetEventList(elb);
                 t->Draw(var+">>h2");
                 
                 // otherwise normalize histos
                 if (norm)
                 {
                         h1.Scale(1./h1.GetEntries());
                         h2.Scale(1./h2.GetEntries());
                 }
                 h1.Scale(Wsig); h2.Scale(Wbkg);
                 
                 // if we want no normalization fill a summed histogram
                 if (!norm) {hs.Add(&h1);hs.Add(&h2);}
                 
                 // set maximum of both histos to display full y-range
                 double maxi = 1.05*max(h1.GetMaximum(), h2.GetMaximum());
                 if (!norm) maxi = hs.GetMaximum()*1.05;
                 h1.SetMaximum(maxi);
                 h2.SetMaximum(maxi);
                 
                 if (!norm) hs.DrawCopy();
                 h1.DrawCopy(norm?"":"same");
                 h2.DrawCopy("same");
                 
                 // draw the cut with direction arrow if requested
                 if (drawcut) 
                 {
                         ln.DrawLine(cutval,0,cutval,maxi*0.9);
                         arr.DrawArrow(cutval,maxi*0.9,cutval+0.1*(xmax-xmin)*cutdir,maxi*0.9,0.01,"|>");
                 }
         }
         
         c1->cd();
 }

bool FCFindLimits	(	TString	var,
		double &	min,
		double &	max
	)

Definition at line 209 of file findcuts.C.

References elpre, h, i, t, and x.

Referenced by FCDrawVariable(), and FCQaVar().

 {
         // consider signal and background after precut
         t->SetEventList(elpre);
         // draw into a dummy histogram; this creates value vectors in the tree
         int nsel = t->Draw(var+">>hhh","","goff");
         t->SetEventList(0);
         
         if (nsel>500000) nsel=500000;
         TH1F *h = (TH1F*)gDirectory->Get("hhh");
         // if no entry left return
         if (h->GetEntries()==0) return false;
         
         // now determine the limits; therefore loop through TTree::GetV1() 
         // being a vector with the values of var
         min=100000., max=-100000.;
         double *x=t->GetV1();
         for (int i=0;i<nsel;++i)
         {
                 if (x[i]>max) max=x[i];
                 if (x[i]<min) min=x[i];
         }
         // increase the upper limit by a tiny value to include this largest value
         // in the uppermost bin, which has a high bin edge '< max', not '<= max'
         max+=1e-6;
         return true;
 }

void FCInit ( int numvars )

Definition at line 349 of file findcuts.C.

References c1, cx, cy, elb, elpre, els, and t.

Referenced by findcuts().

 {
         gStyle->SetTitleFontSize(0.08);
     gStyle->SetPadTopMargin(0.08);
         gStyle->SetOptStat(0);
         
         t->SetEventList(0);
         t->SetBranchStatus("*",1);
 
         // configure the TCanvas division and sizes of histos
         double plotdim = 250;
         if (numvars>11) plotdim=200;
         int cx = numvars+1;
         if (numvars>4)  cx=(numvars+2)/2;
         if (numvars>11) cx=5;
         if (numvars>14) cx=6;
         if (numvars>17) cx=7;
         
         int cy = (numvars+1)/cx;
         if ((numvars+1)%cx) cy+=1;
         
         // refetch canvas or create a new one
         c1=(TCanvas*)gROOT->FindObject("c1");   
         if (c1==0) c1 = new TCanvas("c1","c1",5,5,1000,650);
         c1->Clear();
         c1->Divide(cx, cy);
         c1->SetWindowSize(cx*plotdim, cy*plotdim*1.1);
         
         // refetch event lists if existing or create new ones
         els = (TEventList*)gROOT->FindObject("els");
         if (els==0) els = new TEventList("els");
         
         elb = (TEventList*)gROOT->FindObject("elb");
         if (elb==0) elb = new TEventList("elb");
 
         elpre = (TEventList*)gROOT->FindObject("elpre");
         if (elpre==0) elpre = new TEventList("elpre");
 }

void FCPrepareTree	(	TString	bnames,
		TString	ctlvar
	)

Definition at line 112 of file findcuts.C.

References bnam, FCSplitString(), i, name, printf(), t, and TString.

Referenced by findcuts().

 {
         // Separate names of branches to be considered
         StrVec vecpatt, vecali;
         FCSplitString(bnames," ", vecpatt);
         
         // Enable/disable branches
         t->SetBranchStatus("*",0);
         bnam.clear();  // this array will contain the names of variables (branches and aliases) to be checked
         
         // ******** BRANCHES
         
         // Enable/disable branches according to list of patterns in vecpatt
         for (uint i=0;i<vecpatt.size();++i)
                 if (vecpatt[i].BeginsWith("!")) t->SetBranchStatus(TString(vecpatt[i](1,1000)),0);
                 else t->SetBranchStatus(vecpatt[i],1);
                 
         // for preparing list disable ctlvar
         if (ctlvar!="") t->SetBranchStatus(ctlvar,0);   
         
         // Create a vector with all variables to be considered
         TObjArray* blist = t->GetListOfBranches();
         for (int i=0; i<=blist->GetLast(); ++i)
         {
                 TString name = ((TBranch*)blist->UncheckedAt(i))->GetName();
                 // Add branch name to list, excludes the control variable 
                 if (t->GetBranchStatus(name)>0) bnam.push_back(name); 
         }
                 
         // for sure enable branch with ctlvar
         if (ctlvar!="") t->SetBranchStatus(ctlvar,1);   
         
         // ******** ALIASES
         
         // now we still need to check the list of aliases and enable the branches connected 
         StrIntMap alimap;  // this serves as a marker map for the active aliases
         TList *alilist = t->GetListOfAliases();
         
         // Are there aliases defined?
         if (alilist)
         {       
                 // loop through all patterns
                 for (uint i=0;i<vecpatt.size();++i)
                 {
                         // first create the regular exp to identify the aliases matching the patterns
                         TString patt = vecpatt[i];
                         patt.ReplaceAll("*",".*");
                         int toggle = 1; // switch the alias on (1) or off (0) ?
                         if (patt.BeginsWith("!")) {toggle = 0; patt=patt(1,1000);}
                         TRegexp regali(patt);
                         
                         // save the alias status (on or off) in the map
                         // therefore loop through alias list and toggle all those to either on or off, which match the pattern
                         for (int j=0; j<=alilist->GetLast(); ++j)
                         {
                                 TString aliname = ((TNamed*)alilist->At(j))->GetName();
                                 if (aliname(regali)!="") alimap[aliname] = toggle;
                         }
                 }
                 
                 // loop through map to save the enabled aliases in the bnam vector
                 StrIntMapIt it = alimap.begin();
                 while (it != alimap.end()) 
                 {
                         if (it->second) { vecali.push_back(it->first); bnam.push_back(it->first); }
                         it++;
                 }
                 
                 // now we finally need to enable all branches connected to one of the enabled aliases
                 for (uint i=0;i<vecali.size();++i)
                 {
                         // get the formula to the enabled alias
                         TString fml = t->GetAlias(vecali[i]);
                         //cout << "Alias "<<vecali[i]<<" -> "<<fml<<endl;
                         TRegexp regvar("[A-Za-z_][A-Za-z_0-9]*");
                         while (fml(regvar)!="")
                         {
                                 TString sbr = fml(regvar);
                                 if (t->GetBranch(sbr)) {t->SetBranchStatus(sbr,1);}
                                 fml.Replace(fml.Index(sbr),sbr.Length(),"");
                         }
                 }
         }
         
         // ******** REPORT
         
         // Print out all branches found
         cout<<"**** Name pattern: \""<<bnames<<"\""<<endl;
         cout<<"**** "<<bnam.size()<<" branches/aliases found: \n";
         for (uint i=0;i<bnam.size();++i) 
                 printf("%s  ",bnam[i].Data());
 
         cout <<endl;
 }

void FCPrintInfo ( )

Definition at line 68 of file findcuts.C.

Referenced by findcuts().

 {
         cout <<"\nfindcuts: Tool to find good selection for signal vs bkg. Ranks variables according to maximum (S*ws)^2/(S*ws+B*wb) on optionally normalised distributions.\n";
         cout <<"          Double click on a pad applies corresponding cut and re-iterates.\n\n";
         cout <<"USAGE:\nfindcuts(TTree *t, TString ctlvar, TString sigcut, TString bnames, TString precut, int numvars, double qaopt, double ws, double wb, int norm, int bins)\n"; 
         cout <<"  t       : tree containing signal and background; if given as only argument, list of available branches are printed.\n";
         cout <<"  ctlvar  : control variable to check signal quality (e.g. invariant mass); displayed in upper left pad; excluded from branch list.\n";
         cout <<"  sigcut  : cut to isolate signal; background is selected with !(sigcut).\n";
         cout <<"  bnames  : blank separated list of branch/alias variable names to be considered; can make use of name* / *name /*name*; !(*)name(*) excludes variables (default = \"*\").\n";
         cout <<"  precut  : cut to be applied before variable ranking is done (default = \"\").\n";
         cout <<"  qaopt   : optimisation mode; qaopt=0 : best S/sqrt(S+B); -1<qaopt<0 : best signal eff for fixed bkg reduction; 0<qaopt<1 : best bkg reduction for fixed signal eff. (default = 0).\n";
         cout <<"  numvars : number of variables to be displayed (default = 9).\n";
         cout <<"  ws, wb  : weight factors for signal & background (default = 1.0); are applied to spectra and calculations after optional normalisation.\n";
         cout <<"  norm    : normalisation of distributions before ranking (-> pdf) (default = 1); double click in control variable pad toggles 'norm' and reruns.\n";
         cout <<"  bins    : number of histogram bins, precision depending on binning (default = 500).\n\n";
         cout <<"Example 0: Print all available branches/aliases in tree:\n";
         cout <<"  root [0] findcuts(mytree)\n\n"; 
         cout <<"Example 1: Analyse all branches in TTree 'mytree', control variable 'mass', and signal cut 'issignal':\n";
         cout <<"  root [0] findcuts(mytree, \"mass\", \"issignal\")\n\n";
         cout <<"Example 2: Analyse all branches expect 'evcnt', 'run', those ending with 'px', 'py', or 'pz' or containing the word 'fit', after precut 'goodtracks>6&&goodphotons>2'\n";
         cout <<"  root [0] findcuts(mytree, \"mass\", \"issignal\", \"* !evcnt !run !*p[xyz] !*fit*\", \"goodtracks>6&&goodphotons>2\")\n\n";
         return;
 }

double FCQaVar	(	TString	var,
		double	qaopt = `0`
	)

Definition at line 247 of file findcuts.C.

References BINS, bkgeff, elb, els, fabs(), FCFindLimits(), FCSfc(), h1, h2, i, normalize, qamap, sigeff, t, TString, Wbkg, Wsig, xmax, and xmin.

Referenced by findcuts().

 {
         double xmin, xmax;
         
         if (FCFindLimits(var, xmin, xmax))
         {
                 TH1F h1("h1","",BINS,xmin,xmax);
                 TH1F h2("h2","",BINS,xmin,xmax);
 
                 // Fill histograms for var
                 
                 // signal histo
                 t->SetEventList(els);
                 t->Draw(var+">>h1","","goff");
                 
                 // background histogram
                 t->SetEventList(elb);
                 t->Draw(var+">>h2","","goff");
                 
                 // normalize histograms (if not significance optimisation, do anyways!)
                 if (normalize || fabs(qaopt)>0.0001)
                 {
                         h1.Scale(1./h1.GetEntries());
                         h2.Scale(1./h2.GetEntries());
                 }
                 
                 // multiply histos with weights (only for significance optimisation)
                 if (fabs(qaopt)<0.0001) {h1.Scale(Wsig); h2.Scale(Wbkg);}
                 double h1ent = h1.Integral(), h2ent = h2.Integral();
                 
                 // some vars to store the best efficiencies, cuts, partial sums, etc
                 double rsums=0., lsums=0., rsumb=0., lsumb=0.;
                 double bestSr = 0., bestSl=0., bestCutl=0., bestCutr=0., besteffsl=0., besteffbl=0., besteffsr=0., besteffbr=0.;
                 
                 // integrate from both sides
                 for (int i=1;i<=BINS;++i)
                 {
                         // sum the bin contents from left
                         lsums+=h1.GetBinContent(i);   // signal
                         lsumb+=h2.GetBinContent(i);   // background
                         
                         // sum the bin contents from right
                         rsums+=h1.GetBinContent(BINS-i+1);  // signal
                         rsumb+=h2.GetBinContent(BINS-i+1);  // background
                         
                         double Sl=0, Sr=0;
                         // compute the significances; weights are already taken into account by scaling the histograms
                         if (abs(qaopt)<1e-8) 
                         {
                                 Sl = FCSfc(lsums,lsumb);
                                 Sr = FCSfc(rsums,rsumb);
                         }
                         // compute the signal efficiency for a certain background suppression
                         else if (qaopt<0)
                         {
                                 if (lsumb<=(1.+qaopt)) Sl = lsums;
                                 if (rsumb<=(1.+qaopt)) Sr = rsums;
                         }
                         // compute the signal efficiency for a certain background suppression
                         else
                         {
                                 if (lsums>=qaopt) Sl = 1.-lsumb;
                                 if (rsumb>=qaopt) Sr = 1.-rsumb;
                         }
                         
                         // find the best ones
                         if (Sl>bestSl) {bestSl=Sl; besteffsl=lsums/h1ent; besteffbl=lsumb/h2ent; bestCutl=h1.GetBinLowEdge(i+1);}
                         if (Sr>bestSr) {bestSr=Sr; besteffsr=rsums/h1ent; besteffbr=rsumb/h2ent; bestCutr=h1.GetBinLowEdge(BINS-i+1);}
                 }
                 
                 // create cut strings
                 TString lcut = TString::Format("%s<%f",var.Data(),bestCutl);
                 TString rcut = TString::Format("%s>%f",var.Data(),bestCutr);
                 
                 // and store in map with 1/(QA value) as key (since map is sorted in increasing order)
                 double invqal = 1./bestSl, invqar = 1./bestSr;
                 
                 // if a cut with the same qa exists, insignificantly increase 
                 // current values slightly for being able to store them; otherwise map entries are replaced
                 int maxsearch = 100;
                 while (qamap.find(invqal)!=qamap.end() && --maxsearch>0) invqal = 1./(1./invqal+0.0001);
                 
                 maxsearch = 100;
                 while (qamap.find(invqar)!=qamap.end() && --maxsearch>0) invqar = 1./(1./invqar+0.0001);
                 
                 // store the qa -> cut combination for left and right cut
                 qamap[invqal] = lcut;
                 qamap[invqar] = rcut;
                 
                 // store the corresponding efficiency values
                 sigeff[lcut] = besteffsl; 
                 bkgeff[lcut] = besteffbl;
                 sigeff[rcut] = besteffsr; 
                 bkgeff[rcut] = besteffbr;
         }
 
         return 0;
 }

double FCSfc	(	double	S,
		double	B
	)

Definition at line 239 of file findcuts.C.

References sqrt().

Referenced by FCQaVar(), and findcuts().

 {
         return S/sqrt(S+B);
 }

int FCSplitString	(	TString	s,
		TString	delim,
		StrVec &	toks
	)

Definition at line 94 of file findcuts.C.

References i, and TString.

Referenced by FCPrepareTree().

 {
         toks.clear();
         TObjArray *tok = s.Tokenize(delim);
         int N = tok->GetEntries();      
         for (int i=0;i<N;++i) 
         {
                 TString token = (((TObjString*)tok->At(i))->String()).Strip(TString::kBoth);
                 toks.push_back(token);
         }
         return toks.size();
 }

int findcuts	(	TTree *	theTree = `0`,
		TString	ctlvar = `""`,
		TString	sigcut = `""`,
		TString	bnames = `""`,
		TString	precut = `""`,
		int	numvars = `9`,
		double	qaopt = `0.`,
		double	ws = `1.`,
		double	wb = `1.`,
		int	norm = `1`,
		int	bins = `500`
	)

Definition at line 499 of file findcuts.C.

References BINS, bkgeff, bnam, c1, cnt, comcurrent, comtemplate, cut, effb, effs, elb, elpre, els, fabs(), FCDrawVariable(), FCInit(), FCPrepareTree(), FCPrintInfo(), FCQaVar(), FCSfc(), i, max(), N0b, N0s, Nb, normalize, Ns, printf(), qamap, sigeff, t, TString, Wbkg, and Wsig.

 {
         // if no argument given, print usage information
         if (theTree==0) {FCPrintInfo();return;} 
         t=theTree;
         
         // check for reasonable qaopt parameter -1<qaopt<1
         if (qaopt<-1) qaopt=-1;
         if (qaopt>1)  qaopt=1;
         cout <<"\nOptimising for ";
         if (qaopt<0) cout <<"best signal effciency for given background reduction of "<<fabs(qaopt)<<endl;
         if (qaopt>0) cout <<"best background reduction for given signal efficiency of "<<qaopt<<endl;
         else cout <<"best significance S/sqrt(S+B)"<<endl;
         
         // set defaults for some parameters
         if (numvars==0) numvars = 9;
         
         // print some status information about the TTree
         TString tname = t->GetName();
         cout <<"\nTTree           : "<<tname<<"\nEvents total    : "<<t->GetEntriesFast()<<"\nBranches total  : "<<t->GetNbranches()<<"\nAliases total   : ";
         cout<< (t->GetListOfAliases() ? t->GetListOfAliases()->GetSize() : 0) <<endl;
         
         // the current command string; will be printed at the end
         comcurrent = TString::Format("findcuts(%s,\"%s\",\"%s\",\"%s\",\"%s\",%d,%.3f,%.3f,%.3f,%%d,%d)",
                                                                                  tname.Data(),ctlvar.Data(), sigcut.Data(),bnames.Data(),precut.Data(),numvars, qaopt, ws, wb, bins);
         
         // The template string for pasting the new cut in
         comtemplate = TString::Format("findcuts(%s,\"%s\",\"%s\",\"%s\",\"%s%sNEWCUT\",%d,%.3f,%.3f,%.3f,%d,%d)",
                                                                 tname.Data(), ctlvar.Data(), sigcut.Data(),bnames.Data(),precut.Data(),precut==""?"":"&&",numvars, qaopt, ws, wb, norm, bins);
         
         // catch some special cases
         // if no sigcut is provided, only the list of branches is printed out
         bool printonly = (sigcut=="");
         
         // empty string for branch names defaults to all branches -> '*'
         if (bnames=="") bnames="*";
         if (precut=="") precut="1";
         Wsig = ws; Wbkg = wb;
         
         // rank with normalized spectra; else use real event numbers
         normalize = norm;
         BINS = bins;
         
         // prepare signal and background cut
         TString bkgcut = "!("+sigcut+")";
         TString sfull  = precut+"&&"+sigcut;
         TString bfull  = precut+"&&"+bkgcut;
                 
         // prepare event lists for signal and background
         if (!printonly)
         {
                 // refetch canvas and eventlists
                 FCInit(numvars);
                 
                 t->Draw(">>els",sfull);
                 t->Draw(">>elb",bfull);
                 t->Draw(">>elpre",precut);
                 
                 Ns  = els->GetN();
                 Nb  = elb->GetN();
                 N0s = t->GetEntries(sigcut); 
                 N0b = t->GetEntries(bkgcut); 
         
                 cout <<"Events selected : "<<elpre->GetN()<<endl;
         }
         cout <<endl;
         
         // select variables and fill StrVec bnam with selected branch names
         // ctlvar will be excluded from the ranking
         FCPrepareTree(bnames,ctlvar);           
         
         if (printonly) return;
         
         if (t->GetBranch(ctlvar)==0) {cout <<"**** Unknown control variable: "<<ctlvar<<endl<<endl; return;}
         
         // Draw the control variable in 1st pad
         if (ctlvar!="") FCDrawVariable(ctlvar,1,0,"ctlvar: "+ctlvar);
         
         // if at least one event signal and background is left do ranking
         if (Nb>0 && Ns>0) 
         {
                 // prepare the map for ranking
                 qamap.clear();
                 
                 // print progress bar and prepare some variables
                 cout <<"+-----------+-----------+\n";
                 double ticks = 0; 
                 
                 // **** do the actual ranking procedure by analysing all variables
                 for (uint i=0;i<bnam.size();++i) 
                 {
                         if (bnam[i]!=ctlvar) FCQaVar(bnam[i], qaopt);
                         // print progress bar
                         while ((double)i/(double)bnam.size()>ticks/25.) { cout <<"#"<<flush; ticks++; }
                 }
                 cout <<endl<<endl;
                 
                 // print out and display the results. qamap is sorted by its key being the qa value itself!
                 DblStrMapIt it = qamap.begin();
                 int cnt=0;
                 
                 while (it != qamap.end() && cnt++<numvars) 
                 {
                         double  qa = it->first;
                         TString cut = it->second;
                 
                         if (sigeff[cut]<=bkgeff[cut]) continue;
                 
                         printf("(%2d) %-25s : qa = %6.3f   (eff_s = %5.3f, eff_b = %5.3f)\n", cnt, cut.Data(), 1./qa, sigeff[cut], bkgeff[cut]);
                         
                         // isolate the variable name from the cut and plot
                         TString var = cut(0,max(cut.Index(">"),cut.Index("<")));                                
                         FCDrawVariable(var,cnt+1,normalize,cut);
                         
                         ++it;
                 }
         }
         else cout <<"\n**** No "<<(Nb==0?"background":"signal")<<" events left. No ranking possible."<<endl;
         
         t->SetBranchStatus("*",1);
         t->SetEventList(0);
         
         // the current command string; will be printed at the end
         cout <<"\nCommand:\n"<<TString::Format(comcurrent.Data(),normalize)<<endl;
         
         // print efficiency values
         double effs = (double)Ns/N0s, effb = (double)Nb/N0b;
         double qa  = normalize ? FCSfc(ws*effs, wb*effb) : FCSfc(ws*Ns, wb*Nb);
         double pur = normalize ? ws*effs/(ws*effs+wb*effb) : ws*Ns/(ws*Ns+wb*Nb);
         printf("\nEFF_S = %5.3f (%d/%d)   EFF_B = %5.3f (%d/%d)   PUR = %5.3f   S/sqrt(S+B) %s= %6.3f\n",effs, Ns, N0s, effb, Nb, N0b, pur, normalize?"norm'd ":"", qa);
         
         c1->cd(); c1->Update();
   return 0;
 }