Installing Open MPI 1.6.5 (Ubuntu 12.04, 13.04, fedora)

The OpenMPI logo. A box filled with useful tricks!
Open MPI (http://www.open-mpi.org )  is one of the most liberating tools out there. In a world where time is of the essence and most computers these days have more than one core, why let them sit around idle doing nothing, put them to good use!!

It stands for Message Parsing Interface (MPI) and allows you to do a whole manner of parallelised computing applications. For an in-depth idea of what it can do visit the website but a few key things that are useful to know are:

  1. Collective functions: Main example is the MPI Reduce function, which allows you to perform simple operations such as a summation using the family of processors to do so.
  2. Point-to-Point communication: Most common example I can think of is the use of a head node splitting a dataset into smaller memory chunks based on the number of sub-processors that are available, where each will then compute the same task in parallel. This operations protocol is usually called a master-to-slave process.
  3. Communicators: These connect all the MPI processes into groups, called process groups. An example is the world communicator which contains attributes such as the size (number of processors) and rank (ordered topology) of the group. They also allow for the manipulation of process groups.

How to Install
Open MPI is relatively simple to install, I should point out I have not tried this on a Fedora machine but as long as you have the same libraries/dependencies then it should be the same procedure. There are two methods; the first (and not my preferred method!) is to simply find the latest version from the on-line repositories as of today the latest version was 1.6.5. To do this you can

sudo apt-get update
sudo apt-get install -y autotools-dev g++ build-essential openmpi1.5-bin openmpi1.5-doc libopenmpi1.5-dev
# remove any obsolete libraries
sudo apt-get autoremove

That should get you what you need, if you are on Fedora simply “sudo yum search openmpi” that should bring you up something similar to openmpi, openmpi-devel. The next method and my preferred way as you get the most recent update version(currently 1.6.5) is to take it directly from the website. The below script worked for me on Ubuntu 13.04 tested on 25/07/2013.

The script below will install Open MPI in your /usr/local area, this can be modified by changing the parameter installDIR in the script to the desired location. After install the libraries are placed in $installDIR/lib/openmpi and you can now begin playing with Open MPI. One thing to note is that I apply ldconfig to the /usr/local/lib, this is a much better method than setting paths explicitly. To do this you need to modify your ld.so.conf.d directory to make it look in the /usr/local/xxx area if it doesn’t already. Now with Ubuntu you may already have this linked up so check if your machine has a file called “/etc/ld.so.conf.d/libc.conf” and a path explicitly showing “/usr/local/lib” and that should be all. If you do then you can ignore this step, else you can add the path using the following:

sudo echo "/usr/local/lib" >> /etc/ld.so.conf.d/local.conf
sudo ldconfig -v

I prefer this method as you do not have to keep adding things to your LD_LIBRARY_PATH all the time which is not really recommended, see http://www.xahlee.info/UnixResource_dir/_/ldpath.html for a couple of examples for the case against setting this path. I should have mentioned this in previous blogs too!!

# Matthew M Reid. install open mpi shell script

# install destination
installDIR="/usr/local"
# First get necessary dependencies
sudo apt-get update
sudo apt-get install -y gfortran autotools-dev g++ build-essential autoconf automake 
# remove any obsolete libraries
sudo apt-get autoremove

# Build using maximum number of physical cores
n=`cat /proc/cpuinfo | grep "cpu cores" | uniq | awk '{print $NF}'`

# grab the necessary files
wget http://www.open-mpi.org/software/ompi/v1.6/downloads/openmpi-1.6.5.tar.gz
tar xzvf openmpi-1.6.5.tar.gz
cd openmpi-1.6.5

echo "Beginning the installation..."
./configure --prefix=$installDIR
make -j $n
make install
# with environment set do ldconfig
sudo ldconfig
echo
echo "...done."

So finally to test the installation here is a simple example that just prints out some information from each processor.

#include <iostream>
#include <mpi.h>

int main(int argc, char *argv[])
{
    int numprocessors, rank, namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocessors);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Get_processor_name(processor_name, &namelen);

    if ( rank == 0 )
    {
        std::cout << "Processor name: " << processor_name << "\n";
	std::cout << "master (" << rank << "/" << numprocessors << ")\n";
    } else {
        std::cout << "slave  (" << rank << "/" << numprocessors << ")\n";
   }
   MPI_Finalize();
   return 0;
}

When compiling C++ with Open MPI you can use the executable wrapper, mpic++, which makes compiling much easier. On execution of your script you need to call mpirun where you can specify the number of processors via the -np flag. The output is as follows running from my local machine. Dont worry about the ordering in which these are spit back at you.

matt@matt-W250ENQ-W270ENQ:$ mpic++ -W -Wall test.cpp -o test.o
matt@matt-W250ENQ-W270ENQ:$ mpirun -np 4 test.o
Slave   (3/4)
Processor name; matt-W!112332-NZ10
Master  (0/4)
Slave   (1/4)
Slave   (2/4)

The thing I would like to come onto will be Boost.MPI. This is a very nice interface to the MPI framework and also allows the use of prerequisite parallel graph libraries, which have been well developed and implemented. So the next blog will be about installing Boost which alone has a vast amount to offer, followed by some examples of the Boost.MPI framework in action. The thing that is really clever about this is the Boost.Serialization architecture which allows you to send more complex data structures such as user defined classes via the MPI framework, so you can make “almost” anything parallel.

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