// --------------------------------------------------------------------------------------
//
// ******** LOAD with '.L findcuts.C+' to use as function in CINT ********
//
// 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.
// Double click on a pad applies corresponding cut and re-iterates.
//
// USAGE: findcuts(TTree *t, TString ctlvar, TString sigcut, TString bnames, TString precut, int numvars, double ws, double wb, int norm)
//
// 	  t       : tree containing signal and background; if given as only argument, list of available branches are printed.
// 	  ctlvar  : control variable to check signal quality (e.g. invariant mass); displayed in upper left pad; excluded from branch list.
// 	  sigcut  : cut to isolate signal; background is selected with !(sigcut).
// 	  bnames  : blank separated list of branch variable names to be considered; can make use of name* / *name /*name*; !(*)name(*) excludes variables (default = '*').
// 	  precut  : cut to be applied before variable ranking is done (default = '').
// 	  numvars : number of variables to be displayed (default = 9).
// 	  ws, wb  : weight factors for signal & background (default = 1.0); are applied to spectra and calculations after optional normalisation.
// 	  norm    : normalisation of distributions before ranking (-> pdf) (default = 1.0); double click in control variable pad toggles and reruns.
//    bins    : number of histogram bins, precision depending on binning (default = 500)
//
// --------------------------------------------------------------------------------------

#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"

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;

TCanvas *c1;                    // canvas to display variables
TEventList *els, *elb, *elpre;  // event lists for signal/background/all
TTree *t;                       // the TTree

int N0s, N0b;                   // initial signal/background count
int Ns, Nb;                     // signal/background count after preselection
double Wsig, Wbkg;              // signal/background weights for computing S^2/(S+B)
TString comtemplate, comcurrent;// template command string for the TExec; and current command string
bool normalize;                 // flag for possible normalisation
int BINS;                       // histogram bins

StrVec bnam;                    // vector containing the branch names to be analysed
DblStrMap qamap;                // maps QA value -> cut
StrDblMap sigeff, bkgeff;       // maps cut -> signal/background efficiency

// -------------------------------------------
// Routine to print out the usage of this tool
void FCPrintInfo()
{
	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;
}

// -------------------------------------------
// Splits a TString into vector of strings; separation character contained in delim 
int FCSplitString(TString s, TString delim, StrVec &toks)
{
	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();
}

// -------------------------------------------
// Prepares the list of branches under consideration
// bnames contains branch list with jokers like *name or !*name* et
// branch ctlvar will be enabled for plotting, but excluded from the ranking list
// also takes into account aliases of TTree
void FCPrepareTree(TString bnames, TString ctlvar)
{
	// 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;
}

// -------------------------------------------
// Determines the minimum and maximum value for variable named var for current selecion (for histograms)
bool FCFindLimits(TString var, double &min, double &max)
{
	// 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;
}

// -------------------------------------------
// Compute significance
double FCSfc(double S, double B)
{
	return S/sqrt(S+B);
}

// -------------------------------------------
// The QA routine. It finds the best value ws*ws*S*S/(ws*S +wb*B) for cuts from left and right
// It stores the found cut strings in the qamap
double FCQaVar(TString var, double qaopt=0)
{
	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;
}

// -------------------------------------------
// Initialize style and find TCanvas if existing
// and set the event lists
void FCInit(int numvars)
{
	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");
}

// -------------------------------------------
// Displays variable in pad #numpad in standard canvas
void FCDrawVariable(TString var, int numpad, int norm=1, TString cut="")
{
	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();
}

// -------------------------------------------
// The main function
// -------------------------------------------
void 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)
{
	// 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();
}