Installing nvidia CUDA toolkit – Linux (Ubuntu 13.10)

NVidia

CUDA is one of Nvidia’s greatest gifts for the everyday programmer who dreams of a parallelised world. With most computer/laptops now taking aboard at least one dedicated graphics processor, the prospect of parallel tools for all is becoming a reality. Soon enough everyone no matter what your sex, will be able to multi-task.

I write this as I recently had to upgrade my laptop. I wanted to reinstate my GeForce 640M Nvidia graphics card installing with the latest drivers (331.49) and software for CUDA (v5.5). This turns out is not so trivial as CUDA-5.5 will not run with the latest gcc compilers (gcc4.6+) at present, so you can’t use C++11 :(… However, I came across the two following posts which provided some novel solutions to this problem and now everything works great! I hope they may help someone else out there. Read these both first, just so you are aware of what need be done. The first is merely the process to install the drivers, the second provides more detail into what you should do to make it work on the latest gcc versions although it is a bit of a hack but does the job. Read them both!

How-to-manually-install-latest-nvidia-drivers-on-ubuntu

Incompatibility with GCC 4.7+

You can check the install was ok by running a simple example and/or typing the following command.

nvidia-smi

I will get around to posting more about CUDA in particular, the use of Thrust which is a STL like interface enabling you to devise neat and familiar C++ looking parallel code for a large chunk of parallel applications . It is seriously cool!

For a neat intro to parallel programming see these talks, they know their stuff:-

https://www.udacity.com/course/viewer#!/c-cs344/l-55120467/m-65830481

C++ neural networks applied to stock markets

stock market moves

I have quite an interest in random processes and exactly how we can go about looking for patterns in seemingly random phenomena. One interesting topic I wanted to bring up today was the use of neural networks to predict future stock market trends. What I will present is some free C++ code you can install and then try for yourselves whereby you can employ your own neural networks (NN) on real market data using some wrappers to ROOT (http://root.cern.ch/drupal/)! The strategy is based on buying a certain stock, in this case it is actually trading with a benchmark index known as the SNP500 (http://us.spindices.com/indices/equity/sp-500). The code is an extension to a wonderful package called Hudson written by Alberto Giannetti, he offers a GPL licensed end of day back testing framework using C++ stl, boost, ta-libs and gnu-gsl libraries ( original code base is here http://code.google.com/p/hudson/ ) where you can load in various formats of csv file (Yahoo, Google etc). This extension allows one to perform a multivariate analysis (MVA) which can decipher linear and non-linear trends between input variables of a given dataset and then test its performance. Please bear in mind that this was written during my PhD studies and was solely for fun one weekend. There are clear areas that can be improved in terms of the C++ I just don’t have the time right now, feel free to email me or leave a comment if you like for ways to aid the project. My extension can be obtained from the following location

# Download the repo
git co https://mattreid9956@bitbucket.org/mattreid9956/hudson-mva.git
cd hudson-mva
./build.sh
cd example/MVA_SNP500
# open the run.sh script and try executing the lines one by one to see if they all work, let me know if you have trouble.

This week there is an accompanying pdf file where you can find the detail of a study employed on SNP500 data obtained from Yahoo Finance.

TMVA-StockMarket.pdf

Anyway check out the pdf file if you find it remotely interesting give me a shout. It’s not “complete”, but close enough and I wanted to get it out there and see if anyone is interested in collaborating (I am trying to finish my PhD as well as several other things and don’t like leaving things on the back burner for too long, so here it is). Anyone up to date with particle physics techniques will immediately recognise the format here, anyone who is not will hopefully gain an insight into how some of these techniques can be applied in a toy model trading situation. A network such as this may not be valid indefinitely, its effectiveness would need continuous evaluation, since like I say in the paper, a market is not static and there will be periods of non-normal activity (government intervention for instance). What would be interesting is developing methods of when to turn on a trend following algorithm and when not to! Anyway, if you have any thoughts let me know!

Python Google Trends API

google

EDIT: FYI PEOPLE —> I RETIRED THE REPO AT BIT BUCKET, not to panic though, dreyco676 has a version for python 3 and it is working well as of 21/10/2014. Please see https://github.com/dreyco676/pytrends.

Hi once again, slight detour but thought I should share this. I found the original code on-line and with a few tweaks, managed to get it to do what I wanted and so here it is for anyone interested in this sort of thing. In a world where large datasets are becoming ever more available to the average Joe, Google are doing their bit by allowing you to see what historic rates of search terms have occurred over a given time period, essentially allowing one to look back at what people have been thinking about. You can try this for yourself by following this link to http://www.google.co.uk/trends homepage. As an example lets say I was curious to see how often people search for the term fruit, I would get the following display.

googletrend

(ARGH, since this blog is free I cannot embed the link!!) There are plenty of ways one can then analyses these forms of data whether it be sentiment indicators such as the stock market, movie hits based on search results, when are most babies born and other such seasonal traffic patterns just as some examples.

Anyway the code is simple python script and follows a simple example so check it out. It does require you to login to your Google account so that it can cache the cookies so if you don’t have one you can always get each “.csv” file you need directly from the Google Trends website posted above. EDIT 20/01/2014: You will not be able to login if you have Google additional security running, this is because you get a redirect java session that will wait for a password that is sent to you by mobile, thus the script will never know what that is and is not written in a way to accept it as input. Don’t turn your security off to use this script, thats just stupid, instead just try/make a different gmail account, sorted!!

# Download the git repo.
git clone https://github.com/dreyco676/pytrends.git
cd pytrends
./example.py

EDIT: –> As mentioned above please use the pytrends version on github. I will not be supporting mine. There is an example called example.py, run this and you sould download a search for “pizza”!

This simple example will go and grab a load of trend data provided in the python list and store each in a .csv file. You can also check out the repository directly at https://bitbucket.org/mattreid9956/google-trend-api. The formatting is such that you are returned the end of week date for the whole week and the trend value over that period, this is supposed to make life easier should one run an analysis later. One could easily modify this script to get the desired formatting, note that period in which you search will change the granularity of the time window. For instance searches for 3 months will return daily results, where as searches over a year will return the accumulated results over a given week. This is a little annoying an I don’t see why Google won’t allow daily results by default, maybe time to ask them! Watch this space…

RooFit – Blind Analysis and Scripts

It is not just High Energy Particle Physics (HEP) that find this kind of analysis useful but of course coming from this background, I will take this view point. The world is full of influence and thus bias. Just as you walk down the high street and observe a world of advertisements all there to influence a decision you may or may not make that day. The fact you have seen it means you make a decision (whether subliminal or not).

So why the hell does anyone wish to carry out a blind analysis? When I first heard of it I thought it was a waste of time, however now I am beginning to appreciate it for what it is worth. A nice paper by P. Harrison from the University Of Warwick summarises the importance in HEP and can be found at http://www.ippp.dur.ac.uk/old/Workshops/02/statistics/proceedings/harrison.pdf. In essence, a blind analysis is a measurement which is performed whilst removing ones self from observing the result until the analysis is practically complete. Blind analyses are the optimal way to reduce or eliminate experimenter’s bias, which is the unintended biasing of a result in a particular direction.
In HEP we wish to get the most accurate systematic error estimates for our measurements, this is difficult or impossible if one of the errors is observer bias. This has been shown historically for measurements such as the speed of light where in the 19th century, one measurement after another recorded results that were far from the “true” value but followed very closely with the previous observation. I mean, if everyone else is measuring it to be one thing and you come out with something that is several standard deviations away then you begin to question your result, and maybe your ability… In any standard HEP analysis one would normally apply some selection criteria (“”cuts””) to discard uninteresting data or events, in order to enhance your sensitivity to the signal you’re trying to measure. If you have direct access to the effect that these cuts have on you signal yield, you may be inclined to choose cut values that affect the answer in some subtle way. Well, enter blind analyses where you shield yourself from the final result, thus removing the chance that you introduce a bias into your analysis.
The main purpose of this blog was to highlight the tools available in RooFit (part of the ROOT data analysis framework) and provide some scripts on how you can blind signal yields. There are two scripts; one is a simple example of how to blind yourself from the yield when fitting a Gaussian function, the second is a more complex example based on a “real” application. Now I wont discuss all the individual aspects of the code outlined here, but if there are questions feel free to comment! So to get the code either use wget as per below or go to the bottom of this post and copy->paste.

# Download and then execute the first Gaussian script.
wget https://dl.dropbox.com/u/88131281/blindGaussian.C
wget https://dl.dropbox.com/u/88131281/fit2CrystalBall.C
root -l blindGaussian.C

The first script here is called blindGaussian.C and should be copied with the same file name. It can be run from the command line, the -l flag means ROOT starts in quiet mode so it doesn’t display its promotional information which is really annoying. For those who are not familiar with ROOT check out my previous blog . This example uses the RooUnblindPrecision class, part of the RooBlindTools. This class transforms the variable you wish to blind so you never see the end result even though it will still go into the fit function as a parameter. This is useful for telling you if your fitting pdf is stable and whether it converges (most often tested by running multiple toy simulations, or “toys” as they are more affectionately known).

Gaussian Fit with blinded yields.
First plot showing the fit to a simple Gaussian function, 10,000 events are sampled from the Gaussian function and RooFit performs an extended unbinned maximum likelihood fit to the data points. The yield of course in known and we would expect 10,000 events. Try it and see what you get before and after unblinding. The isBlind flag sets the option to hide the output, it is clear that in this very simple example one generates 10,000 events thus one expects a fit to the data to be of that order. What you should notice is that with blinding the yield is 9.9742e-01, which is clearly mad. If you remove the flag, thus the unblind fit is performed and you return 1.0000e+04, with the statistical errors and everything, wonderful…

Next the more interesting plots. The next piece of code I wanted to show is of course FAKE but a more realistic example aimed at physicists. These types of plots come about a lot in conferences and it is not often I have found really good examples online of how to go about making them. I hope those who are interested will find a working example with which will save them time and effort. Here I construct a Probability Density Function (PDF) made up of two signal like crystal-ball functions, followed by a second order Chebychev polynomial to model background.

Blinded out both Bs and Bd mass regions.
First plot shows the blinded regions, both the B_d and B_s are removed from the plot and the number of events, or yields, are transformed hiding the true values.

Full mass range with individual pdf plotted.
Second plot shows the full unblinded fit with the combined pdf (blue) and the 3 composite pdfs that make went into making it.

Full mass range with fit residuals plotted along the bottom.
Final plot is as the previous but the residuals comparing the distance between the constructed pdf and the data points is compared. Generally a decent fit in this case and should give a \chi^2/ndof \approx 1, moreover a better figure of merit is the probability which is calculated at the very end of the script using root TMath package. A probability, P(\chi^2,ndof) of 50% is a perfect scenario and means that you have neither over or under-estimated the errors. If you are interested in the mathematics of the $\chi^2$ test let me know I will consider doing another blog on that and how to optimise a binning scheme, but you can find some information on Wikipedia, http://en.wikipedia.org/wiki/Chi-squared_distribution.

// Blinding the signal yield of a Gaussian fit.
void blindGaussian()
{
  // Set a flag to blind or unblind, be cautious of this in practice you don't want 
  // to unblind yourself by mistake...
  Bool_t isBlind(kTRUE);
  //Bool_t isBlind(kFALSE);
  
  // Expected yield, this parameter sets the max number of events to be generated in
  // the Gaussian, hence when we fit we "expect" to see a number of this magnitude 
  // returned in our final result.
  double expectedYield(10000);

  // First declare the regular variables
  RooRealVar x("x","x", -10.0, 10.0) ;
  RooRealVar sigma("sigma","sigma", 1.0, 0.1, 2.0) ;
  RooRealVar mean("mean","mean", 0.0) ;
  mean.setConstant();

  // RooCategory used for blind/unblind switching.
  TString blind("blind"), unblind("unblind");
  RooCategory blindCat("blindCat","blind state Category");
  blindCat.defineType(unblind, 0);
  blindCat.defineType(blind, 1);
  
  if(isBlind)
      blindCat.setLabel(blind);
  else 
      blindCat.setLabel(unblind);

  // Generate a dataset with a true m of 0.70
  RooGaussian gaussian("gaussian","gaussian", x, mean, sigma) ;
  RooDataSet* data = gaussian.generate(x, expectedYield) ;
  
  // Define the parameter containing the signal yields from the fit. Always allow floatation
  // below 0, in practice yields can be negative. I allow a range of 2 * the maximum expected, 
  // in general should stop the errors calculations reaching limits.
  RooRealVar nsig("nsig", "nsig", expectedYield/2.0, -expectedYield * 2.0, expectedYield * 2.0);

  // Create an unblinding transformation (a function object deriving from RooAbsReal).
  // This object will return the unblind value corresponding to the blind value in nsig.
  // In this example we use the 'precision blinding' technique as found in
  // RooBlindTools, which is a direct translation of the BlindTools package.
  // Please consult the BlindTools documentation for details on the blinding technique itself
  RooUnblindPrecision m_unblind_nsig("m_unblind_nsig","nsig (unblind)", "BlindString", nsig.getVal(), 100, nsig, blindCat, 0);
  
  // Build a PDF feeding it the unblinded value of m instead of m itself
  RooAbsPdf* extended = new RooExtendPdf("extened","extended", gaussian, m_unblind_nsig); 
  
  // Fit data with gfit (using a blind deltam)
  RooFitResult* fitResult = extended->fitTo(*data, RooFit::Extended(kTRUE), RooFit::NumCPU(2), RooFit::Save(kTRUE), RooFit::Minos(kTRUE)) ;

  // m_unblind will not reveal its contents easily
  m_unblind_nsig.Print() ;
  fitResult->Print("V");
  TCanvas* c1 = new TCanvas("c1", "c1", 900, 500);
  RooPlot* frame = x.frame( RooFit::Title("Best Plot Ever"));
  data->plotOn(frame, RooFit::Binning(40), RooFit::DataError( RooFit::RooAbsData::SumW2 ) );
  extended->plotOn(frame, RooFit::LineColor(4), RooFit::Normalization( 1.0, RooFit::RooAbsReal::Relative));
  frame->Draw();
  frame->SetMinimum(1.e-5); 
  c1->Draw();

  delete extended; extended=0;
}
// Faked data for a blinded fit to two crystal ball functions
#include <math.h>
#include <string>
#include <sstream>
	
using namespace RooFit;

//-------------------------------------------------------------------------------
//-------------------------- Configure parameters -------------------------------
// Master blinding tool, if false everything is unblinded
Bool_t isBlind(kTRUE);
// Define blinding Bs or Bd regions
Bool_t blindBs(kTRUE);
Bool_t blindBd(kTRUE);
Bool_t isResidual(kFALSE);


// Set fileName and paths
TString massVariable("B_M_JpsiConstr");
TString decayMode("B #rightarrow J/#psi K_{S}^{0}");
Int_t bins(40);
// Mass Ranges
Double_t B_mass_lowest(5167.0), B_mass_highest(5478.0), 
	 Bd_mass_min(5251.0), Bd_mass_max(5309.0),
	 Bs_mass_min(5330.0), Bs_mass_max(5393.0);
// Estimates of expected yields for signal and background
Double_t E_Bd_nsig(700.0), E_Bs_nsig(400.0), E_background(5000.0);
//-------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
//
// FUNCTION DEFINITION:
//
//-------------------------------------------------------------------------------

// Evaluate the residual:
double residual( double datum, double pdf )
{
	double chi2 = 0.;

	if ( pdf > 0 )
		chi2 += 2. * ( pdf - datum );

	if ( datum > 0 && pdf > 0 )
		chi2 += 2. * datum * log( datum / pdf );

	return ( ( datum >= pdf ) ? sqrt( chi2 ) : -sqrt( chi2 ) );
}

// Make the RooHist of the residual.
TH1D* residualHist( const RooHist* rhist, const RooCurve* curve )
{
	double r = 0.2;
	double sr = 1. / r;

	// Grab info from the histogram.
	int     n = rhist->GetN();
	double* x = rhist->GetX();
	double* y = rhist->GetY();

	// Create residual histogram.
	double xMin = x[ 0     ];
	double xMax = x[ n - 1 ];
	TH1D* residuals_temp = new TH1D( "r", "", n, xMin, xMax );

	double datum = 0.;
	double pdf   = 0.;

	// Fill the histogram.
	if ( curve )
		for ( int bin = 0; bin < n; bin++ )
		{
			datum = y[ bin ];
			pdf   = curve->Eval( x[ bin ] );
			residuals_temp->SetBinContent( bin + 1, residual( datum, pdf ) );
			residuals_temp->SetBinError  ( bin + 1, 1. );
		}

	residuals_temp->SetMinimum    ( -5.   );
	residuals_temp->SetMaximum    (  5.   );
	residuals_temp->SetStats      ( false );
	residuals_temp->SetMarkerStyle( 8     );
	residuals_temp->SetMarkerSize ( .8    );

	TAxis* xAxis = residuals_temp->GetXaxis();
	xAxis->SetTickLength ( sr * xAxis->GetTickLength()  );
	xAxis->SetLabelSize  ( sr * xAxis->GetLabelSize()   );
	xAxis->SetTitleSize  ( sr * xAxis->GetTitleSize()   );
	xAxis->SetLabelOffset( sr * xAxis->GetLabelOffset() );

	TAxis* yAxis = residuals_temp->GetYaxis();
	yAxis->SetNdivisions ( 504                          );
	yAxis->SetLabelSize  ( sr * yAxis->GetLabelSize()   );

	return residuals_temp;
}

// To draw lines in the residual histograms:
TLine* midLine;
TLine* uppLine;
TLine* lowLine;

void Lines( double xMin, double xMax )
{
	midLine = new TLine( xMin,  0., xMax,  0. );
	uppLine = new TLine( xMin,  2., xMax,  2. );
	lowLine = new TLine( xMin, -2., xMax, -2. );

	uppLine->SetLineColor( kRed );
	lowLine->SetLineColor( kRed );

	midLine->Draw( "same" );
	uppLine->Draw( "same" );
	lowLine->Draw( "same" );
}

void Lines( const TH1D* resid )
{
	double xMin = resid->GetXaxis()->GetXmin();
	double xMax = resid->GetXaxis()->GetXmax();

	midLine = new TLine( xMin,  0., xMax,  0. );
	uppLine = new TLine( xMin,  2., xMax,  2. );
	lowLine = new TLine( xMin, -2., xMax, -2. );

	uppLine->SetLineColor( kRed );
	lowLine->SetLineColor( kRed );

	midLine->Draw( "same" );
	uppLine->Draw( "same" );
	lowLine->Draw( "same" );
}

TString convertToString(Double_t number)
{
	std::ostringstream ss;
	ss <<  number;
	//return std::to_string(number);
	return TString(ss.str());
}

Double_t convertToDouble( const std::string& s )
{
	std::istringstream i(s);
	double x;
	if (!(i >> x))
		return 0;
	return x;
} 

//-------------------------------------------------------------------------------
//
// Main fitting sequence
//
//-------------------------------------------------------------------------------
void fit2CrystalBall()
{
	gROOT->SetStyle("Plain");

	if(isBlind && (!blindBs && !blindBd))
	{
		std::cout << "WARNING:"<< std::endl;
		std::cout << "    isBlind is true, but you have a region you wish to unblind."<< std::endl;
		std::cout << "    for clarity blindBd=" << blindBd << " and blindBs=" << blindBs << std::endl;
		return;
	}
	
	RooRealVar  *MASS   = new RooRealVar(massVariable, "M_{"+decayMode+"}", B_mass_lowest, B_mass_highest, "MeV/c^{2}");

	Double_t dmass(5.), Bs_mean(5366.77), Bd_mean(5279.58);
	RooRealVar *cbmean_core = new RooRealVar("cbmean_core", "cbmean core", Bd_mean, Bd_mean - dmass, Bd_mean + dmass) ; cbmean_core->setConstant(kTRUE);
	RooRealVar *cbsigma_core = new RooRealVar("cbsigma_core", "cbsigma core" , 6, 1, 10.0) ;
	RooRealVar *n_core = new RooRealVar("n_core","", 1.15304);
	RooRealVar *alpha_core = new RooRealVar("alpha_core","", 2.32123);
	RooAbsPdf *cball_core = new RooCBShape("cball_core", "crystal ball core", *MASS, *cbmean_core, *cbsigma_core, *alpha_core, *n_core);

	RooRealVar *cbmean_tail = new RooRealVar("cbmean_tail", "cbmean tail", Bd_mean, Bd_mean - dmass, Bd_mean + dmass); cbmean_tail->setConstant(kTRUE);
	RooRealVar *cbsigma_tail = new RooRealVar("cbsigma_tail", "cbsigma tail" , 5, 1, 10.0) ;
	RooRealVar *n_tail = new RooRealVar("n_tail","", 0.396915);
	RooRealVar *alpha_tail = new RooRealVar("alpha_tail","", -3.11431);
	RooAbsPdf *cball_tail = new RooCBShape("cball_tail", "crystal ball tail", *MASS, *cbmean_core, *cbsigma_tail, *alpha_tail, *n_tail);

	RooRealVar *fmass_core = new RooRealVar("fmass_core", "Bd mass core fraction", 0.9,0.,1.0);
	RooAddPdf  *Bmass_signal = new RooAddPdf("Bmass_signal","Bd mass core+tail gauss", RooArgList(*cball_core,*cball_tail), *fmass_core);

	RooRealVar *Bscbmean_core  = new RooRealVar("Bscbmean_core", "Bscbmean core", Bs_mean, Bs_mean - dmass, Bs_mean + dmass); Bscbmean_core->setConstant(kTRUE);
	RooRealVar *Bscbsigma_core = new RooRealVar("Bscbsigma_core", "Bscbsigma core", 5, 1, 10.0) ;
	RooRealVar *Bsn_core       = new RooRealVar("Bsn_core","", 1.15304);
	RooRealVar *Bsalpha_core   = new RooRealVar("Bsalpha_core","", 2.32123);
	RooAbsPdf  *Bscball_core   = new RooCBShape("Bscball_core", "Bs crystal ball core", *MASS, *Bscbmean_core, *Bscbsigma_core, *Bsalpha_core, *Bsn_core);

	RooRealVar *Bscbmean_tail  = new RooRealVar("Bscbmean_tail", "Bs cbmean tail", Bs_mean, Bs_mean - dmass, Bs_mean + dmass); Bscbmean_tail->setConstant(kTRUE);
	RooRealVar *Bscbsigma_tail = new RooRealVar("Bscbsigma_tail", "Bs cbsigma tail", 5, 1, 10.0) ;
	RooRealVar *Bsn_tail       = new RooRealVar("Bsn_tail","", 0.396915);
	RooRealVar *Bsalpha_tail   = new RooRealVar("Bsalpha_tail","", -3.11431);
	RooAbsPdf  *Bscball_tail   = new RooCBShape("Bscball_tail", "Bs crystal ball tail", *MASS, *Bscbmean_core, *Bscbsigma_tail, *Bsalpha_tail, *Bsn_tail);

	RooRealVar *Bsfmass_core  = new RooRealVar("Bsfmass_core", "Bd mass core fraction", 0.9,0.,1.0);
	RooAddPdf  *Bsmass_signal = new RooAddPdf("Bsmass_signal","Bd mass core+tail gauss", RooArgList(*Bscball_core,*Bscball_tail),*Bsfmass_core);

	RooRealVar c0("c0","c0 parameter", 0.4, -1., 1.) ;
	RooRealVar c1("c1","c1 parameter", -0.07, -1., 1.) ;
	RooAbsPdf *pdf_back = new RooChebychev("pdf_back","Background", *MASS, RooArgList(c0,c1));

	// Generate a data set from the pdfs	
  	RooDataSet* nsigdata = Bmass_signal->generate( *MASS, E_Bd_nsig) ;
  	RooDataSet* Bsnsigdata = Bsmass_signal->generate( *MASS, E_Bs_nsig) ;
  	RooDataSet* bkgdata = pdf_back->generate( *MASS, E_background);
	RooDataSet *SignalData = new RooDataSet("data","data", RooArgSet(*MASS));
	SignalData->append(*nsigdata); SignalData->append(*Bsnsigdata); SignalData->append(*bkgdata);

	// To stop fit hitting bounds make sure to float yields with large range AND negative
	Double_t Nentries(static_cast<Double_t>SignalData->sumEntries());
	RooRealVar nsig("nsig","B_{d}^{0} Signal Yield", E_Bd_nsig, -2. * Nentries, 2. * Nentries);
	RooRealVar Bsnsig("Bsnsig","B_{s}^{0} Signal Yield", E_Bs_nsig, -2. * Nentries, 2. * Nentries);	
	RooRealVar nbkg("nbkg","Background Yield", E_background, -2. * Nentries, 2. * Nentries);

	/////////////////////////////////////
	///// Blind Procedure
	/////////////////////////////////////	

	// RooCategory used for blind/unblind switching.
	TString blind("blind"), unblind("unblind");
	RooCategory blindCatBd("blindCatBd","Bd blind state Category");
	RooCategory blindCatBs("blindCatBs","Bs blind state Category");
	blindCatBd.defineType(unblind, 0);
	blindCatBs.defineType(unblind, 0);
	blindCatBd.defineType(blind, 1);
	blindCatBs.defineType(blind, 1);
	if(isBlind && blindBd)
		blindCatBd.setLabel(blind);
	else 
		blindCatBd.setLabel(unblind);
	if(isBlind && blindBs)
		blindCatBs.setLabel(blind);
	else 
		blindCatBs.setLabel(unblind);
	
	// Unblinding transformation
	RooUnblindPrecision B_nsig("B_nsig","nsig B0 blind","nsigB0blinded", nsig.getVal(), 1000.,  nsig, blindCatBd);
	RooUnblindPrecision B_Bsnsig("B_Bsnsig","nsig Bs blind","nsigBsblinded", Bsnsig.getVal(), 1000., Bsnsig, blindCatBs);

	RooAbsPdf *ext_pdf_sig   = new RooExtendPdf("ext_pdf_sig","Ext. PDF Signal",*Bmass_signal, B_nsig);
	RooAbsPdf *ext_pdf_Bssig = new RooExtendPdf("ext_pdf_Bssig","Ext. PDF Signal",*Bsmass_signal, B_Bsnsig); 	
	RooAbsPdf *ext_pdf_bkg     = new RooExtendPdf("ext_pdf_bkg","Ext. PDF back",*pdf_back, nbkg);

	// The total extended pdf
	RooAbsPdf *pdf_total = new RooAddPdf("pdf_total","Total PDF", RooArgList(*ext_pdf_sig, *ext_pdf_Bssig, *ext_pdf_bkg));

	
	MASS->setRange("lowersideband", B_mass_lowest, Bd_mass_min);
	MASS->setRange("uppersideband", Bs_mass_max, B_mass_highest);
	MASS->setRange("unblind_signal_Bd", B_mass_lowest, Bs_mass_min);
	MASS->setRange("middleregion", Bd_mass_max, Bs_mass_min);
	MASS->setRange("unblind_signal_Bs", Bd_mass_max, B_mass_highest);
	TCut lowersideband = TCut(massVariable + "<" + convertToString(Bd_mass_min));
	TCut uppersideband = TCut(TString( massVariable + ">" + convertToString(Bs_mass_max) ) );
	TCut middleregion = TCut(TString( massVariable + "<" + convertToString(Bs_mass_min) + "&&" + massVariable + ">" + convertToString(Bd_mass_max) ) );
	TCut unblind_signal_Bd = TCut(TString( massVariable + "<" + convertToString(Bs_mass_min) ) );
	TCut unblind_signal_Bs = TCut(TString( massVariable + ">" + convertToString(Bd_mass_max) ) );

	RooFitResult* fitResult = pdf_total->fitTo(*SignalData, NumCPU(4), Extended(kTRUE), Save(kTRUE), Minos(kFALSE));   
	RooMsgService::instance().setSilentMode(kTRUE);

	// Plot PDF and toy data overlaid
	TCanvas* canvas = new TCanvas("canvas","B_{d}^{0} Mass pdf",1000,700) ;

	if(isResidual) canvas->SetWindowSize(1000, 800);
	Double_t chi2(0.);

	Double_t nd1(0.), nd2(0.), nd3(0.);
	RooPlot *p3 = MASS->frame(Title("Invariant Mass Plot M_{"+decayMode+"}"));
	if(!isBlind){
		SignalData->plotOn(p3, Binning(bins), Name("data_hist"), DataError( RooAbsData::SumW2 ));  
		pdf_total->plotOn(p3, Range("fitRange"), Normalization(1.0,RooAbsReal::RelativeExpected), Name("main_curve"));
		pdf_total->plotOn(p3, Range("fitRange"), Components(*ext_pdf_bkg),LineStyle(kDashed),LineColor(kMagenta),Normalization(1.0,RooAbsReal::RelativeExpected));
		pdf_total->plotOn(p3, Range("fitRange"), Components(*ext_pdf_sig),LineStyle(kDashed),LineColor(kRed),Normalization(1.0,RooAbsReal::RelativeExpected));   
		pdf_total->plotOn(p3, Range("fitRange"), Components(*ext_pdf_Bssig),LineStyle(kDashed),LineColor(kGreen),Normalization(1.0,RooAbsReal::RelativeExpected));   	
	}
	else
	{
		if(blindBd && blindBs)
		{
			std::cout << "Keeping both Bs and Bd regions blind!" << std::endl;
			nd1=(static_cast<Double_t>((SignalData->reduce(lowersideband)->numEntries())))/(static_cast<Double_t>(Nentries)); 
			nd2=(static_cast<Double_t>((SignalData->reduce(middleregion)->numEntries())))/(static_cast<Double_t>(Nentries)); 
			nd3=(static_cast<Double_t>((SignalData->reduce(uppersideband)->numEntries())))/(static_cast<Double_t>(Nentries)); 

			SignalData->plotOn(p3, CutRange("lowersideband"), Name("data_hist1"), Binning(bins) );
			pdf_total->plotOn( p3, Range("lowersideband"), Normalization(nd1,RooAbsReal::Relative), LineWidth(3), LineColor( kBlue), Name("main_curve1"));	
			SignalData->plotOn(p3, CutRange("middleregion"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("middleregion"), Normalization(nd2,RooAbsReal::Relative), LineWidth(3), LineColor( kBlue));	
			SignalData->plotOn( p3, CutRange("uppersideband"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("uppersideband"), Normalization(nd3,RooAbsReal::Relative), LineWidth(3), LineColor( kBlue));
		}
		else if(!blindBd && blindBs)
		{
			std::cout << "Unblinding the Bd mass region, keep Bs blinded!" << std::endl;
			nd1=(static_cast<Double_t>((SignalData->reduce(unblind_signal_Bd)->numEntries())))/(static_cast<Double_t>(Nentries)); 
			nd2=(static_cast<Double_t>((SignalData->reduce(uppersideband)->numEntries())))/(static_cast<Double_t>(Nentries)); 

			SignalData->plotOn(p3, CutRange("unblind_signal_Bd"), Name("data_hist1"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("unblind_signal_Bd"), Normalization(nd1, RooAbsReal::Relative), LineWidth(3), LineColor( kBlue), Name("main_curve1"));
			pdf_total->plotOn(p3, Range("unblind_signal_Bd"), Components(*ext_pdf_sig),LineStyle(kDashed),LineColor(kRed),Normalization(nd1,RooAbsReal::RelativeExpected));   
			SignalData->plotOn( p3, CutRange("uppersideband"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("uppersideband"), Normalization(nd2, RooAbsReal::Relative), LineWidth(3), LineColor( kBlue));		

		}
		else
		{
			std::cout << "Unblinding the Bs mass region, keep Bd blinded!" << std::endl;
			nd1=(static_cast<Double_t>((SignalData->reduce(lowersideband)->numEntries())))/(static_cast<Double_t>(Nentries)); 
			nd2=(static_cast<Double_t>((SignalData->reduce(unblind_signal_Bs)->numEntries())))/(static_cast<Double_t>(Nentries)); 

			SignalData->plotOn(p3, CutRange("lowersideband"), Name("data_hist1"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("lowersideband"), Normalization(nd1,RooAbsReal::Relative), LineWidth(3), LineColor( kBlue), Name("main_curve1"));	
			SignalData->plotOn(p3, CutRange("unblind_signal_Bs"), Binning(bins), DataError( RooAbsData::SumW2 ) );
			pdf_total->plotOn( p3, Range("unblind_signal_Bs"), Normalization(nd2, RooAbsReal::Relative), LineWidth(3), LineColor( kBlue));
			pdf_total->plotOn(p3, Range("unblind_signal_Bs"), Components(*ext_pdf_Bssig),LineStyle(kDashed),LineColor(kGreen),Normalization(nd1,RooAbsReal::RelativeExpected));   
		}
	}
	
	Float_t xmin(0.12), xmax(.37), ymin(.68), ymax(.88);

	if (isResidual && !isBlind)
	{
		p3->GetXaxis()->SetLabelSize(0);
		p3->GetXaxis()->SetTitleOffset(1.15);
		p3->GetYaxis()->SetTitleOffset(1.26);
		RooHist* histogram = p3->getHist("data_hist");
		RooCurve* curve = p3->getCurve("main_curve");
		TH1D* hresidual  = residualHist(histogram,curve);
		canvas->Divide( 1, 2, .1, .1 );
		double r  = .25;
		double sl = 1. / ( 1. - r );
		p3->GetYaxis()->SetLabelSize(sl * p3->GetYaxis()->GetLabelSize());
		//double labelSize = 0.2;
		hresidual->GetXaxis()->SetLabelSize((1./(1.+r)) * hresidual->GetXaxis()->GetLabelSize());
		hresidual->GetYaxis()->SetLabelSize((1./(1.+r)) * hresidual->GetYaxis()->GetLabelSize());
		TPad* padHisto = (TPad*) canvas->cd(1);
		TPad* padResid = (TPad*) canvas->cd(2);
		double small = 0.1;
		padHisto->SetPad( 0., r , 1., 1. );
		padHisto->SetBottomMargin( small );
		padResid->SetPad( 0., 0., 1., r  );
		padResid->SetBottomMargin( 0.3  );
		padResid->SetTopMargin   ( small );
		padHisto->cd();
		TPaveText *box= new TPaveText( xmin, ymin, xmax, ymax, "NDC");
		box->SetFillColor(10);
		box->SetBorderSize(0);
		box->SetTextAlign(12);
		box->SetTextSize(0.04F);
		box->SetFillStyle(1001);
		box->SetFillColor(10);
		Char_t buf[30], buf2[30], buf3[30], bufflhcb[30];
		sprintf( buf,  "B^{0} nsig = %3.2f", nsig.getVal() );
		sprintf( buf2,  "B_{s}^{0} nsig= %3.2f", Bsnsig.getVal() );
		sprintf( buf3,  "nbkg= %3.2f", nbkg.getVal() );
		TText* textbd = box->AddText( buf );
		textbd->SetTextSize(0.035);
		textbd->SetTextAlign(11);
		//textbd->SetTextColor(kRed);
		TText* textbs = box->AddText( buf2 );
		textbs->SetTextSize(0.035);
		textbs->SetTextAlign(11);
		//textbs->SetTextColor(kGreen);
		TText* textbkg = box->AddText( buf3 );
		textbkg->SetTextSize(0.035);
		textbkg->SetTextAlign(11);
		//textbkg->SetTextColor(kMagenta);
		p3->addObject(box); 
		p3->Draw() ;
		//myPave->Draw("SAME");
		padResid->cd();
		hresidual->Draw();
		Lines( hresidual );
		hresidual->Draw( "SAME" );
	}
	else
	{
		p3->GetXaxis()->SetTitleOffset(1.15);
		p3->GetYaxis()->SetTitleOffset(1.26);
		gPad->SetLeftMargin(0.1);
		TPaveText *box= new TPaveText( xmin, ymin, xmax, ymax, "NDC");
		// Or set the exact coordinatres -> TPaveText *box= new TPaveText( 5180., 160., 5260., 200., "br");
		box->SetFillColor(10);
		box->SetBorderSize(0);
		box->SetTextAlign(12);
		box->SetTextSize(0.04F);
		box->SetFillStyle(1001);
		Char_t buf[30], buf2[30], buf3[30];
		sprintf( buf,  "B_{d}^{0} yield = %3.2f", nsig.getVal() );
		sprintf( buf2,  "B_{s}^{0} yield = %3.2f", Bsnsig.getVal() );
		sprintf( buf3,  "bkg yield = %3.2f", nbkg.getVal() );
		TText* textbd = box->AddText( buf );
		textbd->SetTextSize(0.035);
		textbd->SetTextAlign(11);
		//textbd->SetTextColor(kRed);
		TText* textbs = box->AddText( buf2 );
		textbs->SetTextSize(0.035);
		textbs->SetTextAlign(11);
		//textbs->SetTextColor(kGreen);
		TText* textbkg = box->AddText( buf3 );
		textbkg->SetTextSize(0.035);
		textbkg->SetTextAlign(11);
		//textbkg->SetTextColor(kMagenta);
		p3->addObject(box); 
		p3->Draw() ;
	}
	
    p3->SetMinimum(1.e-5);
	canvas->Modify();
	canvas->Draw();
	canvas->SaveAs("MassPlot.eps");
	canvas->SaveAs("MassPlot.png");
	fitResult->Print("V");

    // to get a chi2 value we need to plot the complete pdf, for this as long as we never draw it we are ok.
	RooPlot *p = MASS->frame();
	SignalData->plotOn(p, Binning(bins), Name("data_chi2"), DataError( RooAbsData::SumW2 ));  
	pdf_total->plotOn(p, Range("fitRange"), Normalization(1.0,RooAbsReal::RelativeExpected), Name("curve_chi2"));
	// There is a method within RooPlot called chiSquare(). This returns the reduced chi2 which when multiplied out
    // by the number of bins gives the actual chi2 fit between the . 
    chi2 = p->chiSquare( "curve_chi2", "data_chi2") * bins;
    // The total number of degrees of freedom is defined as the number of bins subtracted by the free floating parameters
    // in the fit. We already know the number of bins as we specified it, so the free params (we also know) but can be
    // obtained using:-
    Int_t floated = static_cast<Int_t>(fitResult->floatParsFinal().getSize());

    std::cout << "Chi2 = " << chi2 << "  (for "<< bins << " bins and " << floated << " floating params)" << std::endl;
    std::cout<<" Prob(chi2,ndof) of above = " << TMath::Prob(chi2, bins-floated) << std::endl;    
    return; 
}