#include <signal.h>

#include <iostream>
#include <sstream>
#include <vector>

#include <cmath>

#include "UtlArguments.hh"
#include "RcdArena.hh"
#include "RcdCount.hh"
#include "RunReader.hh"
#include "SubAccessor.hh"

#include "MpsAnalysis.hh"

#include "MapsTrack.hh"
#include "ToString.h"

MpsAnalysis::MpsAnalysis(unsigned wombat, int runNumber, int bins, int low,
		int high, int trackingSensor) :
	_thresholdScanRun(false), _thresholdRun(false),
			_mpsBeamThresholdScan(false), _mpsBeam(false), _nBt(0),
			_configCount(0) {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	_runNumber = runNumber;
	_bins = bins;
	_lowEdge = low;
	_highEdge = high;
	//Assuming we have N boards, which one is considered for special attention in this code
	_specialSensor = trackingSensor;

	//define (unphysical) sensor Ids for now; they'll get populated by populateSensorIds()
	_sensorIds[0] = 95;
	_sensorIds[1] = 96;
	_sensorIds[2] = 97;
	_sensorIds[3] = 98;
	_sensorIds[4] = 99;
	_pmtDaqIndex = 100;

}

MpsAnalysis::~MpsAnalysis() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	saveToDisk();
	delete _rootFile;
	std::cout << "~MpsAnalysis(): done. Destructing."<< std::endl;
}

void MpsAnalysis::saveToDisk() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	char name[100];
	strcpy(name, "MpsAnalysisTracks_");
	strcat(name, toString(_runNumber).c_str());
	strcat(name, ".root");
	_mtm.setRequiredLeftSensor(_mySensors[8]);
	_mtm.setRequiredRightSensor(_mySensors[6]);

	_mtm.exportToRootFile(name);
//
//	_mtmDispT.setRequiredLeftSensor(_mySensors[8]);
//	_mtmDispT.setRequiredRightSensor(_mySensors[6]);
//	char name2[100];
//	strcpy(name2, "MpsAnalysisTracksDispT_");
//	strcat(name2, toString(_runNumber).c_str());
//	strcat(name2, ".root");
//	_mtmDispT.exportToRootFile(name2);
	writeRootFile();

	_rootFile->Close();
}

void MpsAnalysis::computeEfficiencies() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	double eff(0);
	for (unsigned sensor(0); sensor < 5; sensor++) {
		for (int thresh(0); thresh < _highEdge; thresh++) {
			if (_ConfirmedTracks[sensor][thresh] > 0) {
				eff = static_cast<double>(_ConfirmedTracks[sensor][thresh])
						/ _Tracks[sensor][thresh] * 100;
				_hEffIntegratedAllBts[sensor]->Fill(thresh, eff);
			}
		}
	}
}

void MpsAnalysis::findNoisyPixels(int cut) {
	std::cout << std::cout << __PRETTY_FUNCTION__ << " with cut = "<< cut
			<< std::endl;
	ofstream highpixels;
	highpixels.open("highpixels.txt");
	unsigned numberOfHighPixels(0);

	for (unsigned sensor(0); sensor < 5; sensor++) {
		highpixels << "Sensor pseudo value: \t"<< sensor << std::endl;

		for (unsigned x(0); x < 168; x++) {
			for (unsigned y(0); y < 168; y++) {
				double hits = _hHitPattern[sensor]->GetBinContent(x, y);
				_hThreshEfficacy[sensor]->Fill(hits);
				if (hits > cut) {
					//naughty pixel! you're going dowwwwwnnn....
					highpixels << x << "\t"<< y << "\t"<< hits << "\n";
					numberOfHighPixels++;
				}
			}
		}
		highpixels << std::endl;
	}
	std::cout << "\tFound "<< numberOfHighPixels << " high pixels"<< std::endl;
	highpixels.close();
}

void MpsAnalysis::findHitDistribution() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	for (unsigned sensor(0); sensor < 5; sensor++) {
		for (unsigned x(0); x < 168; x++) {
			for (unsigned y(0); y < 168; y++) {
				double hits = _hHitPattern[sensor]->GetBinContent(x, y);
				_hThreshEfficacy[sensor]->Fill(hits);
			}
		}
	}
}

void MpsAnalysis::makeGlobalTimeStampPlot() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	for (unsigned sensor(0); sensor < 5; sensor++) {
		_hGlobalTimeStamps->Add(_hTimeStampsBySensor[sensor]);
	}
}

bool MpsAnalysis::determineRunType(
		std::vector<const IlcRunStart*> runStartHeader) {

	assert(runStartHeader.size()==1);

	for (unsigned i(0); i<runStartHeader.size(); i++) {
		//vs[i]->print(std::cout) << std::endl;
		_runNumber = runStartHeader[i]->runNumber();
		// Check if it is a MAPS threshold scan run
		_thresholdScanRun=(runStartHeader[i]->runType().type()==IlcRunType::mpsThresholdScan);
		//or one in the beam
		_mpsBeamThresholdScan=(runStartHeader[i]->runType().type()==IlcRunType::mpsBeamThresholdScan);
		//or just a threshold run
		_thresholdRun = (runStartHeader[i]->runType().type()==IlcRunType::mpsThreshold);

		_mpsBeam = (runStartHeader[i]->runType().type()==IlcRunType::mpsBeam);

		_mpsCosmicsThresholdScan = (runStartHeader[i]->runType().type() == IlcRunType::mpsCosmicsThresholdScan);
	}
	
	if(_thresholdScanRun)
		_mpsCosmicsThresholdScan = true;
	//if we are able to do the run, return true...
	if (_thresholdScanRun || _mpsBeamThresholdScan || _thresholdRun || _mpsBeam
			|| _mpsCosmicsThresholdScan)
		return true;
	

	//otherwise we're not set up to handle this
	return false;
}

void MpsAnalysis::populateSensorIds(SubAccessor myAccessor) {
	std::vector<const MpsLocationData<MpsPcb1ConfigurationData>* >
			v(myAccessor.access< MpsLocationData<MpsPcb1ConfigurationData> >());
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	//There are nWrites which preceed a number of reads - we get the
	// sensor Ids from the reads
	unsigned nWrites(0);
	for (unsigned i(0); i<v.size(); i++) {
		if (!v[i]->write()) {
			std::cout<< "\tSensor id for index "<< i << ":\t"
					<< static_cast<unsigned short>(v[i]->location().sensorId()) << std::endl;
			if (v[i]->location().usbDaqMaster()) {
				//this one is the usbDaqMaster and has the PMTs plugged into it.
				//record this here,
				_pmtDaqIndex = i-nWrites;
				std::cout << "\tUsbDaqMaster is on sensor " << v[i]->location().sensorId() << ", index " << i-nWrites << std::endl;
			}

			_sensorIds[i-nWrites] = v[i]->location().sensorId();
		} else {
			nWrites++;
		}
	}
	std::cout << "Found "<< nWrites
			<< " writes in PCB configuration data, with v.size() = "<< v.size()
			<< std::endl;
}

bool MpsAnalysis::record(const RcdRecord& r) {

	//std::cout<< __PRETTY_FUNCTION__<< std::endl;

	//	if (doPrint(r.recordType())) {
	//		r.RcdHeader::print(std::cout, " ") << std::endl;
	//	}

	SubAccessor accessor(r);

	switch (r.recordType()) {

	case RcdHeader::runStart: {

		// Get list of IlcRunStart subrecords
		// There should only be one
		std::vector<const IlcRunStart*> vs(accessor.access<IlcRunStart>());
		determineRunType(vs);
		break;
	}

	case RcdHeader::runEnd: {
		if (_thresholdScanRun) {
			// Get list of IlcRunEnd subrecords
			// There should only be one
			std::vector<const IlcRunEnd*> ve(accessor.access<IlcRunEnd>());
			assert(ve.size()==1);
			for (unsigned i(0); i<ve.size(); i++) {
				//ve[i]->print(std::cout) << std::endl;
			}
		}
		break;
	}

	case RcdHeader::configurationStart: {
		std::vector<const IlcConfigurationStart*>
				vs(accessor.access<IlcConfigurationStart>());
		assert(vs.size()==1);

		// Get a list of PCB configuration data objects
		std::vector<const MpsLocationData<MpsPcb1ConfigurationData>* >
				v(accessor.access< MpsLocationData<MpsPcb1ConfigurationData> >());

		if (_configCount == 0) {
			//get sensor ids
			populateSensorIds(accessor);
			_mySensors = getMapsSensors();
			for (std::map<unsigned, MapsSensor*>::const_iterator it =
					_mySensors.begin(); it != _mySensors.end(); it++) {
				_mtmDispT.addSensor((*it).second);
				_mtm.addSensor((*it).second);
			}
			//make histograms now we know the sensor Ids
			makeHistograms(_sensorIds);
		}
		_configCount++;
		// Get list of IlcConfigurationStart subrecords
		// There should only be one


		unsigned nWrites(0);
		for (unsigned i(0); i<v.size(); i++) {
			//the first one is a write, the five remaining are reads
			//v[i]->print(std::cout) << std::endl;
			if (!v[i]->write()) {
				// Get threshold differences
				_shaperThreshold[i - nWrites]=v[i]->data()->region01ThresholdValue();
				_samplerThreshold[i - nWrites]=v[i]->data()->region23ThresholdValue();
				//				std::cout << "Sensor:\t" << _sensorIds[i-nWrites] << "\n";
				//				std::cout << "\tShaper threshold:\t"
				//						<< _shaperThreshold[i -nWrites] << "\n";
				//				std::cout << "\tSampler threshold:\t"
				//						<< _samplerThreshold[i -nWrites] << "\n";
			} else {
				nWrites++;
			}
		}

		//a mpsBeamThresholdScan includes an efficiency calculation
		if (_mpsBeamThresholdScan) {
			_thresholdScanRun = true;
		}

		std::vector<const MpsLocationData<MpsUsbDaqConfigurationData>* >
				w(accessor.access<MpsLocationData<MpsUsbDaqConfigurationData> >());

		_spillCycleCount = w[0]->data()->spillCycleCount();

		// Initialise hit counter
		_nHits=0;
		_tHits = 0;
		//_nBt = 0;

		break;

	}

	case RcdHeader::configurationEnd: {
		if (_thresholdScanRun) {

			// Get list of IlcConfigurationEnd subrecords
			// There should only be one
			std::vector<const IlcConfigurationEnd*>
					vs(accessor.access<IlcConfigurationEnd>());

			assert(vs.size()==1);

			for (unsigned i(0); i<vs.size(); i++) {
				//vs[i]->print(std::cout) << std::endl;
			}

		}

		break;
	}

	case RcdHeader::bunchTrain: {
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* >
				v(accessor.access<MpsLocationData<MpsSensor1BunchTrainData> >());
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* >
				u(accessor.access<MpsLocationData<MpsUsbDaqBunchTrainData> >());
		handleBunchTrain(v, u);
		break;
	}

	default: {
		break;
	}
	};

	return true;
}

void MpsAnalysis::handleBunchTrain(
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* > v,
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > u) {
	if (_nBt%1000== 0) {
		std::cout << "Processing bunch train: "<< _nBt << std::endl;
	}
	if (_mpsBeam) {
		handleBeamRun(v, u);
	}

	if (_mpsBeamThresholdScan) {
		//is the threshold being varied, or are we at nominal?
		unsigned scannedSensor(0);
		if ((_configCount - 1)%41!= 0) {
			//which sensor is being scanned?
			scannedSensor = static_cast<int>(((_configCount -1)%41- 1)/10);
		}

		handleBeamThresholdScanRun(v, u, scannedSensor);
	}

	if (_thresholdScanRun) {
		handleThresholdScanRun(v, u);
	}

	if (_mpsCosmicsThresholdScan) {
		//is the threshold being varied, or are we at nominal?
		unsigned scannedSensor(0);
		if ((_configCount - 1)%41!= 0) {
			//which sensor is being scanned?
			scannedSensor = static_cast<int>(((_configCount -1)%41- 1)/10);
		}
		handleCosmicThresholdScanRun(v, u, scannedSensor);
	}

	_nBt++;
}

void MpsAnalysis::handleCosmicThresholdScanRun(
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* > v,
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > u,
		unsigned scannedSensor) {

	memset(_timeStampsByShap, 0, 5*8192*sizeof(unsigned));
	memset(_timeStampsBySamp, 0, 5*8192*sizeof(unsigned));
	memset(_xPosBySensor, 0, 5*8192*2*sizeof(double));
	memset(_yPosBySensor, 0, 5*8192*2*sizeof(double));

	std::map<unsigned, double> pos = getSensorPositions();

	//loop over all the sensors - the first is the PMT, so we ignore it
	for (unsigned i(0); i<v.size(); i++) {
		if (!v[i]->write()) {
			std::vector<MpsSensor1Hit> vh(v[i]->data()->hitVector());
			populateTimeStampsAndPositions(i, vh, 100);
		} else {
			std::cout<< "RcdHeader in bunch train is a write? " << std::endl;
		}
	}

	//now we count hits
	for (unsigned t(0); t < 8192; t++) {
		unsigned hitsInOtherSensorsShap(0);
		unsigned hitsInOtherSensorsSamp(0);
		unsigned hitsInScannedSensorShap(0);
		unsigned hitsInScannedSensorSamp(0);

		unsigned makeTrack(0);
		unsigned nTrackHits(0);
		unsigned minHits(1);

		MapsSensor* scanSensor = _mySensors[_sensorIds[scannedSensor]];

		//Now loop over the sensors and count how many hits each sensor has
		for (unsigned sensor(0); sensor < 5; sensor++) {

			_hTimeStampsBySensor[sensor]->Fill(_timeStampsBySensor[sensor][t]);
			if (sensor != scannedSensor && _sensorIds[sensor] != 99) {
				nTrackHits += _timeStampsByShap[sensor][t]
						+ _timeStampsBySamp[sensor][t];
				if (_timeStampsBySamp[sensor][t] == minHits)
					hitsInOtherSensorsSamp++;
				if (_timeStampsByShap[sensor][t] == minHits)
					hitsInOtherSensorsShap++;
				if ((_timeStampsBySamp[sensor][t]
						+ _timeStampsByShap[sensor][t]) == minHits)
					makeTrack++;
			} else if (sensor == scannedSensor) {
				if (_timeStampsBySamp[sensor][t] == minHits)
					hitsInScannedSensorSamp += _timeStampsBySamp[sensor][t];
				if (_timeStampsByShap[sensor][t] == minHits)
					hitsInScannedSensorShap += _timeStampsByShap[sensor][t];
			}
		}

		//how many hits are there in the threshold-scanned sensor
		unsigned nScanned(hitsInScannedSensorSamp + hitsInScannedSensorShap);

		if (makeTrack > 2) {
			//we assume it's a track

			std::map<MapsSensor*, MapsTrack::coord> hits;
			bool suddenDeath(false);

			for (unsigned sensor(0); sensor < 5; sensor++) {
				bool doHitsFor4th(true);
				if (sensor == scannedSensor && nScanned != 1)
					doHitsFor4th = false;

				if (_sensorIds[sensor] != 99) {
					MapsSensor* theSensor = _mySensors[_sensorIds[sensor]];
					MapsTrack::coord c(
							static_cast<unsigned>(_xPosBySensor[sensor][t][0]),
							static_cast<unsigned>(_yPosBySensor[sensor][t][0]));

					if (doHitsFor4th)
						hits[theSensor] = c;
					if (c.first > 168 || c.second > 168)
						suddenDeath = true;

				}
			}
			//if (nTrackHits == 3 || nTrackHits) {
			//std::cout << __PRETTY_FUNCTION__ << __LINE__ << "\n";
			MapsTrack* aTrack = new MapsTrack(t, hits, scanSensor, _shaperThreshold[scannedSensor]);
			std::pair<double, double> error(0.15, 0.15);
			aTrack->setTrackError(error);
			//diagnose(std::cout, *aTrack);
			//std::cout << *aTrack << "\n";
			//std::cout << *track;

			if (nScanned > 1) {
				std::cout << "Fourth sensor has " << nScanned << " hits.\n";
				suddenDeath = true;
			}
			

			if (aTrack->chiSqProb(0) > 0.001 && aTrack->chiSqProb(1) > 0.001) {
				diagnose(std::cout, *aTrack);
				//				std::cout << "\tnScanned:\t" << nScanned << "\n";
				//				std::cout << "\tnScanndDispT:\t" << nScannedDispT << "\n";
				//				if (nScannedDispT > 0) {
				//					std::cout << "*** Displaced Track! ***" << "\n";
				//					diagnose(std::cout, *aDispTrack);
				//				}
			}
			if (!suddenDeath) {
				_mtm.addTrack(aTrack);
			} else {
				//std::cout << "Sudden death called!\n";
				//diagnose(std::cout, *aTrack);
			}
		}
	}

}

void MpsAnalysis::fillPmtTimestampPlot(
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > uPmts) {

	for (unsigned l = 0; l < uPmts.size(); l ++) {
		//extract pointer to first element of PMT datums
		const MpsUsbDaqBunchTrainDatum* pmtDatum(uPmts[l]->data()->data());

		//loop over PMT tags, we'll label them with h
		for (unsigned h = 0; h < uPmts[l]->data()->numberOfTags(); h++) {
			if (pmtDatum[h].channel(0))
				_hPmtTags1->Fill(pmtDatum[h].timeStamp());
			if (pmtDatum[h].channel(1))
				_hPmtTags2->Fill(pmtDatum[h].timeStamp());
			if (pmtDatum[h].channel(0) && pmtDatum[h].channel(1))
				_hPmtCoincidence->Fill(pmtDatum[h].timeStamp());
		}
	}
}

void MpsAnalysis::populateTimeStampsAndPositions(unsigned sensor,
		const std::vector<MpsSensor1Hit>& hitVec, int lowestThresh) {

	for (unsigned j(0); j < hitVec.size(); j++) {
		//_hHitPattern[sensor]->Fill(hitVec[j].pixelX(), hitVec[j].pixelY());
		if (hitVec[j].pixelX() > 83) {
			//sampler
			_timeStampsBySamp[sensor][hitVec[j].timeStamp()]++;
			if (_shaperThreshold[sensor] > lowestThresh)
				_hTimestampsSamplers->Fill(hitVec[j].timeStamp());
		} else {
			//shaper
			_timeStampsByShap[sensor][hitVec[j].timeStamp()]++;
			if (_shaperThreshold[sensor] > lowestThresh)
				_hTimestampsShapers->Fill(hitVec[j].timeStamp());
		}
		//record the sum and sum of squares of these hits at this timestamp
		//std::cout << "t:\t" << vh[j].timeStamp() << ", (x, y):\t" << vh[j].pixelX() << ", \t" << vh[j].pixelY() << std::endl;
		_xPosBySensor[sensor][hitVec[j].timeStamp()][0] += hitVec[j].pixelX();
		_xPosBySensor[sensor][hitVec[j].timeStamp()][1] += hitVec[j].pixelX() * hitVec[j].pixelX();
		_yPosBySensor[sensor][hitVec[j].timeStamp()][0] += hitVec[j].pixelY();
		_yPosBySensor[sensor][hitVec[j].timeStamp()][1] += hitVec[j].pixelY() * hitVec[j].pixelY();
		assert(hitVec[j].timeStamp() < 8192);
	}
	//let's make a cluster
	if (_shaperThreshold[sensor] > lowestThresh) {
		for (unsigned t(0); t < 8192; t++) {
			unsigned nHits(_timeStampsByShap[sensor][t]
					+ _timeStampsBySamp[sensor][t]);
			if (nHits > 1) {
				double meanX = static_cast<double>(_xPosBySensor[sensor][t][0])
						/ nHits;
				double meanY = static_cast<double>(_yPosBySensor[sensor][t][0])
						/ nHits;
				double meanX2 =
						static_cast<double>(_xPosBySensor[sensor][t][1])
								/ nHits;
				double meanY2 =
						static_cast<double>(_yPosBySensor[sensor][t][1])
								/ nHits;
				double sigmaR = sqrt(meanX2 + meanY2 - meanX * meanX - meanY
						*meanY);

				_hNHitsAndDispersion->Fill(sigmaR, nHits);
				if (sigmaR < 2) {
					_hNClustersInT->Fill(t);
				}

			}
		}
	}
}

std::map<unsigned, double> MpsAnalysis::getSensorPositions() {
	std::map<unsigned, double> pos;
	if (_mpsBeamThresholdScan) {
		pos[6] = 54;
		pos[2] = 36;
		pos[7] = 18;
		pos[8] = 0;
	} else if (_mpsCosmicsThresholdScan) {
		pos[2] = 0;
		pos[6] = 12;
		pos[8] = 24;
		pos[7] = 36;
	}
	return pos;
}

std::map<unsigned, MapsSensor*> MpsAnalysis::getMapsSensors() {
	std::cout << __PRETTY_FUNCTION__ << "\n";
	std::map<unsigned, MapsSensor*> sensors;
	std::map<unsigned, double> zpos = getSensorPositions();
	if (_mpsBeamThresholdScan) {
		//initialise sensors
		for (unsigned k(0); k < 4; k++) {
			MapsSensor* sensor = new MapsSensor(_sensorIds[k], zpos[_sensorIds[k]]);
			sensors[_sensorIds[k]] = sensor;
			if (_sensorIds[k] == 2 || _sensorIds[k] == 8)
				sensor->setPhi(3.14159265358979323846);
			std::cout << "\t" << "SensorIds[k]: "<< _sensorIds[k] << ": "
					<< *sensor << "\n";
		}
	} else if (_mpsCosmicsThresholdScan) {
		//initialise sensors
		for (unsigned k(0); k < 4; k++) {
			MapsSensor* sensor = new MapsSensor(_sensorIds[k], zpos[_sensorIds[k]]);
			sensors[_sensorIds[k]] = sensor;
			if (_sensorIds[k] == 6)
				sensor->setPhi(1.5 * 3.14159265358979323846);
			if (_sensorIds[k] == 8)
				sensor->setPhi(3.14159265358979323846);
			if (_sensorIds[k] == 7)
				sensor->setPhi(0.5 * 3.14159265358979323846);

			std::cout << "\t" << "SensorIds[k]: "<< _sensorIds[k] << ": "
					<< *sensor << "\n";
		}
	} else {
		std::cout << __PRETTY_FUNCTION__ << ":\t";
		std::cout
				<< "Cannot initialise MapsSensors as I don't understand the run type!\n";
	}
	return sensors;
}

bool MpsAnalysis::wasSensorRotated(const unsigned sensorId) {
	if (sensorId == 2 || sensorId ==8)
		return true;
	return false;
}

void MpsAnalysis::handleBeamThresholdScanRun(
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* > v,
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > u,
		unsigned scannedSensor) {

	memset(_timeStampsByShap, 0, 5*8192*sizeof(unsigned));
	memset(_timeStampsBySamp, 0, 5*8192*sizeof(unsigned));
	memset(_xPosBySensor, 0, 5*8192*2*sizeof(double));
	memset(_yPosBySensor, 0, 5*8192*2*sizeof(double));

	std::map<unsigned, double> pos = getSensorPositions();

	//loop over all the sensors - the first is the PMT, so we ignore it
	for (unsigned i(0); i<v.size(); i++) {
		if (!v[i]->write()) {
			std::vector<MpsSensor1Hit> vh(v[i]->data()->hitVector());
			populateTimeStampsAndPositions(i, vh, 100);
		} else {
			std::cout<< "RcdHeader in bunch train is a write? " << std::endl;
		}
	}

	//Extract PMT data, for what it's worth
	fillPmtTimestampPlot(u);

	//How many tracks there are in this BT
	unsigned nPossibleTracks(0);

	//how many hits are there in the threshold-scanned sensor, when a track in made at a given BX
	unsigned nScannedVeto(0);
	unsigned nScannedVetoShap(0);
	unsigned nScannedVetoSamp(0);

	//define DispT to mean 4096 BX from the examined BX (forward or backward)
	//so that we can have a fake-noise calculation
	//(4906 BXs should not have correlated behaviour with the current BX)
	unsigned nScannedVetoDispT(0);
	unsigned nScannedVetoShapDispT(0);
	unsigned nScannedVetoSampDispT(0);

	//Legacy code
	unsigned nPossibleTracksNotSpecial(0);

	//now we count hits
	for (unsigned t(0); t < 8192; t++) {
		//hits in sensors at nominal threshold
		unsigned hitsInOtherSensorsShap(0);
		unsigned hitsInOtherSensorsSamp(0);
		//hits in threshold-scanned sensor
		unsigned hitsInScannedSensorShap(0);
		unsigned hitsInScannedSensorSamp(0);

		//how many sensors have hits > minHits so we can make a track
		unsigned makeTrack(0);

		unsigned nTrackHits(0);

		//for the de-correlated hit count
		unsigned hitsInScannedSensorShapDispT(0);
		unsigned hitsInScannedSensorSampDispT(0);

		//for soft photon studies (legacy code)
		unsigned hitsInSpecialSensor(0);
		
		//how many hits are required in each sensor
		unsigned minHits(1);
		
		//initialise displaced time
		unsigned dispT(0);
		
		//extract threshold-scanned sensor
		MapsSensor* scanSensor = _mySensors[_sensorIds[scannedSensor]];

		//Now loop over the sensors and count how many hits each sensor has
		for (unsigned sensor(0); sensor < 5; sensor++) {

			_hTimeStampsBySensor[sensor]->Fill(_timeStampsBySensor[sensor][t]);
			//sensors at constant threshold...
			if (sensor != scannedSensor && _sensorIds[sensor] != 99) {
				//N.B. If there's more than one hit, give up at the moment
				nTrackHits += _timeStampsByShap[sensor][t]
						+ _timeStampsBySamp[sensor][t];
				if (_timeStampsBySamp[sensor][t] == minHits)
					hitsInOtherSensorsSamp++;
				if (_timeStampsByShap[sensor][t] == minHits)
					hitsInOtherSensorsShap++;
				//If == minHits, we can make a track with this sensor at least
				if ((_timeStampsBySamp[sensor][t]
						+ _timeStampsByShap[sensor][t]) == minHits)
					makeTrack++;
			} else if (sensor == scannedSensor) {
				if (_timeStampsBySamp[sensor][t] == minHits)
					hitsInScannedSensorSamp += _timeStampsBySamp[sensor][t];
				if (_timeStampsByShap[sensor][t] == minHits)
					hitsInScannedSensorShap += _timeStampsByShap[sensor][t];
				
				//define some displaced time
				if (t >= 4096)
					dispT = t - 4096;
				else
					dispT = t + 4096;

				if (_timeStampsBySamp[sensor][dispT] == minHits)
					hitsInScannedSensorSampDispT
							+= _timeStampsBySamp[sensor][dispT];
				if (_timeStampsByShap[sensor][dispT] == minHits)
					hitsInScannedSensorShapDispT
							+= _timeStampsByShap[sensor][dispT];

			}
			//hang over from earlier days :-| !! Not required any more I think
			//Legacy code...
			if (sensor == _specialSensor) {
				handleCosmicThresholdScanRun(v, u, scannedSensor);
				hitsInSpecialSensor+= _timeStampsBySamp[sensor][t]
						+ _timeStampsByShap[sensor][t];
			} else {
				nPossibleTracksNotSpecial += _timeStampsBySamp[sensor][t]
						+ _timeStampsByShap[sensor][t];
			}

		}

		//how many hits are there in the threshold-scanned sensor?
		unsigned nScanned(hitsInScannedSensorSamp + hitsInScannedSensorShap);
		//as above, but at the displaced time
		unsigned nScannedDispT(hitsInScannedSensorSampDispT
				+ hitsInScannedSensorShapDispT);
		
		// we want three hits from the three nominal threshold sensors
		if (makeTrack > 2) {
			//we assume it's a track

			std::map<MapsSensor*, MapsTrack::coord> hits;
			std::map<MapsSensor*, MapsTrack::coord> hitsDispT;

			//new histo - fill with threshold, number of hits in sensor four
			if (scannedSensor == _specialSensor) {
				_hSoftPhotonsSpecialSensor->Fill(
						_shaperThreshold[scannedSensor], hitsInSpecialSensor);
			}

			//plot the mean (x,y) position of the hits in the track
			double sumTrackX(0);
			double sumTrackY(0);

			//sum of squares
			double sumTrackX2(0);
			double sumTrackY2(0);
			
			//this will be used as a veto for storing the track
			bool suddenDeath(false);

			for (unsigned sensor(0); sensor < 5; sensor++) {
				//if (sensor != scannedSensor && _sensorIds[sensor] != 99) {
				bool doHitsFor4th(true);
				bool doHitsFor4thDispT(true);
				if (sensor == scannedSensor) {
					//only do 4th sensor's hits if it has hits anyway!
					if (nScanned != 1)
						doHitsFor4th = false;
					if (nScannedDispT != 1)
						doHitsFor4thDispT = false;
				}
				if (_sensorIds[sensor] != 99) {
					MapsSensor* theSensor = _mySensors[_sensorIds[sensor]];
					//All the even id'd sensors were rotated 180degs w.r.t odd id'd sensors
					//But this rotation is now handled by the MapsSensor objects themselves
					unsigned offset((_timeStampsByShap[sensor][t]
							+ _timeStampsBySamp[sensor][t])*168);
					
					unsigned dispToffset((_timeStampsByShap[sensor][dispT]
							+ _timeStampsBySamp[sensor][dispT])*168);

					sumTrackX += _xPosBySensor[sensor][t][0];
					sumTrackY += _yPosBySensor[sensor][t][0];

					sumTrackX2 += _xPosBySensor[sensor][t][0]
							*_xPosBySensor[sensor][t][0];
					sumTrackY2 += _yPosBySensor[sensor][t][0]
							*_yPosBySensor[sensor][t][0];
					//Store the raw hit coordinate information here
					MapsTrack::coord c(
							static_cast<unsigned>(_xPosBySensor[sensor][t][0]),
							static_cast<unsigned>(_yPosBySensor[sensor][t][0]));

					if (sensor == scannedSensor) {
						if (doHitsFor4thDispT) {
							//need to extract displaced spacetime coordinate
							MapsTrack::coord
									cDispT(
											static_cast<unsigned>(_xPosBySensor[sensor][dispT][0]),
											static_cast<unsigned>(_yPosBySensor[sensor][dispT][0]));
							hitsDispT[scanSensor] = cDispT;
						}
					} else {
						hitsDispT[theSensor] = c;
					}
					if (doHitsFor4th)
						hits[theSensor] = c;
					//veto multiple hits
					if (c.first > 168 || c.second > 168)
						suddenDeath = true;

				}
			}
			//if (nTrackHits == 3 || nTrackHits) {
			//std::cout << __PRETTY_FUNCTION__ << __LINE__ << "\n";
			
			//Initialise the tracks
			MapsTrack* aTrack = new MapsTrack(t, hits, scanSensor, _shaperThreshold[scannedSensor]);
			MapsTrack* aDispTrack = new MapsTrack(t, hitsDispT, scanSensor, _shaperThreshold[scannedSensor]);
			//diagnose(std::cout, *aTrack);
			//std::cout << *aTrack << "\n";
			//std::cout << *track;

			if (nScanned > 1) {
				std::cout << "Fourth sensor has " << nScanned << " hits.\n";
				suddenDeath = true;
			}
			if (nScannedDispT > 1) {
				std::cout << "Fourth sensor has (dispT) " << nScannedDispT
						<< " hits.\n";
				suddenDeath = true;
			}
			//If it's a good track, let's hear about it (based on unaligned system
			//make sure rotations and sensor orders are correct, otherwise this will be meaningless
			//You can of course use the MapsTracks realignment tool to reinitialise
			//sensor parameters.
			if (aTrack->chiSqProb(0) > 0.01 && aTrack->chiSqProb(1) > 0.01) {
				diagnose(std::cout, *aTrack);
				//				std::cout << "\tnScanned:\t" << nScanned << "\n";
				//				std::cout << "\tnScanndDispT:\t" << nScannedDispT << "\n";
				//				if (nScannedDispT > 0) {
				//					std::cout << "*** Displaced Track! ***" << "\n";
				//					diagnose(std::cout, *aDispTrack);
				//				}
			}
			//everything's ok, so let's store the track
			if (!suddenDeath) {
				_mtm.addTrack(aTrack);
				_mtmDispT.addTrack(aDispTrack);
			} else {
				//std::cout << "Sudden death called!\n";
				//diagnose(std::cout, *aTrack);
			}
			//}
			//The rest is all legacy/depracted code
			double meanX(sumTrackX / nTrackHits);
			double meanY(sumTrackY / nTrackHits);

			double stdDevX(sqrt(sumTrackX2/nTrackHits -meanX * meanX));
			double stdDevY(sqrt(sumTrackY2/nTrackHits -meanY * meanY));
			_hTracks3VerticesXDispersion->Fill(stdDevX);
			_hTracks3VerticesYDispersion->Fill(stdDevY);
			bool realTrack;

			//apply a cut on tracks that are really due to beam particles
			if (aTrack->theta() < 0.01) {
				realTrack = true;
				_hTracks3Vertices->Fill(nTrackHits);
				//if (nTrackHits < 4)
				nPossibleTracks++;
				if (_shaperThreshold[scannedSensor] > 120) {
					_hTracks3VerticesInT->Fill(t);
					_hTracks3VerticesMeanXY->Fill(sumTrackX / nTrackHits,
							sumTrackY / nTrackHits);
				}
			} else {
				//it's not really a track
				realTrack = false;
			}

			_hTracksXYDispAnd4thSensor->Fill((stdDevX + stdDevY)/2, nScanned);
			_hTracksXYDispAndHitsInTrack->Fill((stdDevX + stdDevY)/2,
					nTrackHits);
			//			std::cout << "NHits: " << nTrackHits << ", mean X: " << sumTrackX
			//					/ nTrackHitstoString(_runNumber).c_str() << ", mean Y: " << sumTrackY / nTrackHits
			//					<< std::endl;

			//does the threshold-scanned sensor confirm a track?
			if (realTrack) {
				if (nScanned > 0) {
					//we're efficient
					//if (nTrackHits < 4)
					nScannedVeto++;
					//where in the Bunch train do we find this?
					if (_shaperThreshold[scannedSensor] > 120) {
						_hTracks4VerticesInT->Fill(t);
						//How many hits confirm this track?
						_hVetoHitCount[scannedSensor]->Fill(nScanned);
						//What's the mean position of the confirming hits?
						double meanX(_xPosBySensor[scannedSensor][t][0]
								/nScanned);
						double meanY(_yPosBySensor[scannedSensor][t][0]
								/nScanned);

						_hVetoHitXY[scannedSensor]->Fill(meanX, meanY, nScanned);
						_hVetoT[scannedSensor]->Fill(t, nScanned);
						//And their std deviations
						_hVetoHitXStdDev[scannedSensor]->Fill(sqrt(_xPosBySensor[scannedSensor][t][1]
								/nScanned - meanX * meanX));
						_hVetoHitYStdDev[scannedSensor]->Fill(sqrt(_yPosBySensor[scannedSensor][t][1]
								/nScanned - meanY * meanY));
					}
				} else {
					// not efficient!
					//assert(_timeStampsByShap[scannedSensor][t]
					//		+ _timeStampsBySamp[scannedSensor][t] == 0);
				}

				if (hitsInScannedSensorShap > 0)
					nScannedVetoShap++;
				if (hitsInScannedSensorSamp > 0)
					nScannedVetoSamp++;
			}
			//now for the decorrelated ones
			if (nScannedDispT > 0)
				nScannedVetoDispT++;
			if (hitsInScannedSensorShapDispT > 0)
				nScannedVetoShapDispT++;
			if (hitsInScannedSensorSampDispT > 0)
				nScannedVetoSampDispT++;

		}

	}

	//Legacy code
	if (nPossibleTracks > 0) {
		//grind through the efficiency calculations
		//eff = number scanned sensor saw / total possible x 100%
		//std::cout << "3+VX tracks:\t" << nPossibleTracks << ", 4VX:\t"
		//		<< nScannedVeto << std::endl;
		_Tracks[scannedSensor][_shaperThreshold[scannedSensor]]
				+= nPossibleTracks;
		_ConfirmedTracks[scannedSensor][_shaperThreshold[scannedSensor]]
				+= nScannedVeto;

		_hEffAndThreshOverall[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVeto)/(nPossibleTracks) * 100);

		_hEffAndThreshSamp[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVetoSamp)/(nPossibleTracks) *100);

		_hEffAndThreshShap[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVetoShap)/(nPossibleTracks) *100);

		_hEffAndThreshOverallDispT[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVetoDispT) /(nPossibleTracks)* 100);

		_hEffAndThreshSampDispT[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVetoSampDispT) /(nPossibleTracks)
						*100);

		_hEffAndThreshShapDispT[scannedSensor]->Fill(
				_shaperThreshold[scannedSensor],
				static_cast<double>(nScannedVetoShapDispT) /(nPossibleTracks)
						*100);

	}

}

void MpsAnalysis::handleBeamRun(
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* > v,
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > u) {

	memset(_timeStampsBySensor, 0, 5*8192*sizeof(unsigned));
	memset(_allTimeStamps, 0, 8192*sizeof(unsigned));
	for (unsigned i(0); i<v.size(); i++) {
		if (!v[i]->write()) {
			std::vector<MpsSensor1Hit> vh(v[i]->data()->hitVector());
			populateTimeStampsAndPositions(i, vh, 0);
		} else {
			std::cout << __PRETTY_FUNCTION__
					<< " : Found a write record in a bunch train!" << std::endl;
		}
	}
	for (unsigned sensor(0); sensor < 5; sensor++) {
		for (unsigned t(0); t < 8192; t++) {
			_hTimeStampsBySensor[sensor]->Fill(_timeStampsBySensor[sensor][t]);
			if (sensor ==0)
				_hSimulTimeStamps->Fill(_allTimeStamps[t]);
		}
	}

	//loop over each timestamp
	unsigned offset(168);
	for (unsigned t(0); t < 8192; t++) {
		unsigned nTrackHits(0);
		//MapsTrack* aTrack = new MapsTrack(t);
		//std::map<unsigned, MapsTrack::coord> hits;
		for (unsigned sensor(0); sensor < 5; sensor++) {
			//define a candidate track by at least 1 hit in each of three 3 sensors
			if (_sensorIds[sensor] != 99) {
				//if there's exclusively ONE hit in this sensor, increase nTrackHits by one
				//if two or more, give up for now :-)
				if (_timeStampsBySensor[sensor][t] == 1) {
					nTrackHits++;
					if (wasSensorRotated(_sensorIds[sensor])) {
						MapsTrack::coord c( static_cast<unsigned>(offset
								- _xPosBySensor[sensor][t][0]),
								static_cast<unsigned>(offset
										- _yPosBySensor[sensor][t][0]));
						//hits[_sensorIds[sensor]] = c;
					} else {
						MapsTrack::coord
								c(
										static_cast<unsigned>(_xPosBySensor[sensor][t][0]),
										static_cast<unsigned>(_yPosBySensor[sensor][t][0]));
						//hits[_sensorIds[sensor]] = c;
					}
				}

			}
		}
		if (nTrackHits > 2) {
			//it's a candidate track
			//save the track
			//aTrack->setHits(hits);
			//_mtUtil.addTrack(aTrack);			

		} else {
			//delete the track!
			//delete aTrack;
		}

	}

}

void MpsAnalysis::handleThresholdScanRun(
		std::vector<const MpsLocationData<MpsSensor1BunchTrainData>* > v,
		std::vector<const MpsLocationData<MpsUsbDaqBunchTrainData>* > u) {
	//Iterate over each sensor labelled by 'i'
	for (unsigned i(0); i<v.size(); i++) {
		if (!v[i]->write()) {

			_hHitsInBtByThreshold[i]->Fill(_shaperThreshold[i], v[i]->data()->numberOfRegionHits(0)+ v[i]->data()->numberOfRegionHits(1));

			for (unsigned theRegion(0); theRegion < 4; theRegion++) {
				//Need to loop over word hits
				unsigned allHitsA(0);
				unsigned allHitsB(0);
				unsigned wordHitsA(0);
				unsigned wordHitsB(0);
				const MpsSensor1BunchTrainDatum* rDatum =v[i]->data()->regionData(theRegion);

				for (unsigned j(0); j < v[i]->data()->numberOfRegionHits(theRegion); j++) {
					for (unsigned c(0); c < 6; c++) {
						if (rDatum[j].channel(c) != 0&& rDatum[j].row() <84) {
							allHitsA++;
						} else if (rDatum[j].channel(c) != 0 && rDatum[j].row()
								> 83) {
							allHitsB++;
						}
					}
					if (rDatum[j].row() < 84)
						wordHitsA++;
					if (rDatum[j].row() > 83)
						wordHitsB++;
					_hGroupsAndRows[i][theRegion]->Fill(rDatum[j].group(),
							rDatum[j].row());
				}
				unsigned thresh(_shaperThreshold[i]);
				if (theRegion == 2|| theRegion == 3) {
					thresh = _samplerThreshold[i];
				}
				_hHits[i][theRegion]->Fill(thresh, allHitsA+ allHitsB);
				_hHitsByRegion[i][theRegion]->Fill(thresh, allHitsA + allHitsB);
				_hWordHits[i][theRegion]->Fill(thresh, wordHitsA + wordHitsB);

			}

			//Loop over hits, count them in an array, 
			//then corroborate with PMT information

			//let the first index of allHits refer to the shapers/samplers, and the third being the PMT
			unsigned allHits[3][8192];
			memset(allHits, 0, 3*8192*sizeof(unsigned));

			if (!_mpsBeamThresholdScan) {
				//as that run type does this itself
				std::vector<MpsSensor1Hit> vh(v[i]->data()->hitVector());

				memset(_timeStampsByShap, 0, 5*8192*sizeof(unsigned));
				memset(_timeStampsBySamp, 0, 5*8192*sizeof(unsigned));
				memset(_xPosBySensor, 0, 5*8192*2*sizeof(double));
				memset(_yPosBySensor, 0, 5*8192*2*sizeof(double));

				populateTimeStampsAndPositions(i, vh, 150);

				//				for (unsigned t(0); t < 8192; t++) {
				//					allHits[0][t] += _timeStampsByShap[i][t];
				//					allHits[1][t] += _timeStampsBySamp[i][t];
				//				}
				//
				//				for (unsigned l = 0; l < u.size(); l ++) {
				//					//extract pointer to first element of PMT datums
				//					const MpsUsbDaqBunchTrainDatum * pmtDatum(u[l]->data()->data());
				//					//loop over PMT tags, we'll label them with h
				//					for (unsigned h = 0; h < u[l]->data()->numberOfTags(); h++) {
				//						if (pmtDatum[h].channel(0))
				//							allHits[2][pmtDatum[h].timeStamp()]++;
				//						if (pmtDatum[h].channel(1))
				//							allHits[2][pmtDatum[h].timeStamp()]++;
				//					}
				//				}
				fillPmtTimestampPlot(u);
			}
		} else {
			std::cout<< "RcdHeader in bunch train is a write? " << std::endl;
		}
	}
}

void MpsAnalysis::makeHistograms(int sensorIds[]) {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	std::string rootFileTitle(("MpsAnalysis_" + toString(_runNumber)+ ".root").c_str());
	_rootFile = new TFile(rootFileTitle.c_str(),"RECREATE");
	_hNHitsPerConfig
			= new TH2I("hNHitsPerConfig", "Number of hits by threshold", _bins, _lowEdge, _highEdge, 30, 0, 30);
	std::cout << "\tInitialising histograms..."<< std::endl;
	std::cout << "\tbins: "<< _bins << ", lowEdge: "<< _lowEdge
			<< ", _highEdge: "<< _highEdge << std::endl;

	_hHitsInRow
			= new TH1F("hHitsInRow","Hits in a row over all BTs, for one thresh;Pixel",42,0,42);

	_hTimestampsSamplers
			= new TH1F("hTimestampsSamplers", "hTimestampsSamplers;Timestamp", 8192, 0, 8192);
	_hTimestampsShapers
			= new TH1F("hTimestampsShapers", "hTimestampsShapers;Timestamp", 8192, 0, 8192);
	globals.push_back(_hTimestampsSamplers);
	globals.push_back(_hTimestampsShapers);

	_hGlobalTimeStamps = new TH1F("hGlobalTimeStamps", "hGlobalTimeStamps", 40, 0, 40);
	globals.push_back(_hGlobalTimeStamps);

	_hSimulTimeStamps = new TH1F("hSimulTimeStamps", "hSimulTimeStamps", 40, 0, 40);
	globals.push_back(_hSimulTimeStamps);

	memset(_timeStampsBySensor, 0, 5*8192*sizeof(unsigned));
	memset(_allTimeStamps, 0, 8192*sizeof(unsigned));
	memset(_Tracks, 0, 5*2010*sizeof(unsigned));
	memset(_ConfirmedTracks, 0, 5*2010*sizeof(unsigned));

	_hSoftPhotonsSpecialSensor
			= new TH2F("hSoftPhotonsSpecialSensor", "Hits in non-tracking front sensor", _bins, _lowEdge, _highEdge, 100, 0, 100);
	efficiencies.push_back(_hSoftPhotonsSpecialSensor);

	_hTracks3Vertices
			= new TH1F("hTracks3Vertices", "Number of hits in tracks with 3 vertices", _bins, _lowEdge, _highEdge);
	efficiencies.push_back(_hTracks3Vertices);

	_hTracks3VerticesInT
			= new TH1F("hTracks3VerticesInT", "Track (3VX) formation in time", 8192/128, 0, 8192);
	efficiencies.push_back(_hTracks3VerticesInT);
	_hTracks4VerticesInT
			= new TH1F("hTracks4VerticesInT", "Track (4VX) formation in time", 8192/128, 0, 8192);
	efficiencies.push_back(_hTracks4VerticesInT);
	_hTracks3VerticesMeanXY
			= new TH2F("hTracks3VerticesMeanXY", "hTracks3VerticesMeanXY",168, 0, 168, 168, 0, 168);
	efficiencies.push_back(_hTracks3VerticesMeanXY);

	_hTracks3VerticesXDispersion
			= new TH1F("hTracks3VerticesXDispersion", "hTracks3VerticesXDispersion", 200, 0, 200);
	_hTracks3VerticesYDispersion
			= new TH1F("hTracks3VerticesYDispersion", "hTracks3VerticesYDispersion", 200, 0, 200);
	efficiencies.push_back(_hTracks3VerticesXDispersion);
	efficiencies.push_back(_hTracks3VerticesYDispersion);

	_hTracksXYDispAnd4thSensor
			= new TH2F("hTracksXYDispAnd4thSensor", "hTracksXYDispAnd4thSensor", 200, 0, 200, 5, 0, 5);
	_hTracksXYDispAndHitsInTrack
			= new TH2F("hTracksXYDispAndHitsInTrack", "hTracksXYDispAndHitsInTrack", 200, 0, 200, 10, 3, 13);
	efficiencies.push_back(_hTracksXYDispAnd4thSensor);
	efficiencies.push_back(_hTracksXYDispAndHitsInTrack);

	//define a hit pattern for each of 5 possible sensors
	for (unsigned s(0); s < 5; s++) {
		_hEffIntegratedAllBts[s]
				= new TH1F(("hEffIntegratedAllBts" + toString(sensorIds[s])).c_str(), ("Efficiency integrated over all BTs, sensor " + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		efficiencies.push_back(_hEffIntegratedAllBts[s]);

		_hTracks3VerticesInXY[s] = new TH2F(("hTracks3VerticesInXY" + toString(sensorIds[s])).c_str(), ("hTracks3VerticesInXY" + toString(sensorIds[s])).c_str(),168, 0, 168, 168, 0, 168);
		efficiencies.push_back(_hTracks3VerticesInXY[s]);

		_hTracks4VerticesInXY[s] = new TH2F(("hTracks4VerticesInXY" + toString(sensorIds[s])).c_str(), ("hTracks4VerticesInXY" + toString(sensorIds[s])).c_str(),168, 0, 168, 168, 0, 168);
		efficiencies.push_back(_hTracks4VerticesInXY[s]);

		_hVetoHitCount[s] = new TH1F(("hVetoHitCount" + toString(sensorIds[s])).c_str(),("hVetoHitCount" + toString(sensorIds[s])).c_str(), 100, 0, 100);
		efficiencies.push_back(_hVetoHitCount[s]);

		_hVetoHitXY[s] = new TH2F(("hVetoHitXY" + toString(sensorIds[s])).c_str(), ("hVetoHitXY" + toString(sensorIds[s])).c_str(), 168, 0, 168, 168, 0, 168);
		_hVetoHitXStdDev[s] = new TH1F(("hVetoHitXStdDev" + toString(sensorIds[s])).c_str(), ("hVetoHitXStdDev" + toString(sensorIds[s])).c_str(), 100, 0, 20);
		_hVetoHitYStdDev[s] = new TH1F(("hVetoHitYStdDev" + toString(sensorIds[s])).c_str(), ("hVetoHitYStdDev" + toString(sensorIds[s])).c_str(), 100, 0, 20);
		_hVetoT[s] = new TH1F(("hVetoT" + toString(sensorIds[s])).c_str(), ("hVetoT" + toString(sensorIds[s])).c_str(), 8192, 0, 8192);

		efficiencies.push_back(_hVetoHitXY[s]);
		efficiencies.push_back(_hVetoHitXStdDev[s]);
		efficiencies.push_back(_hVetoHitYStdDev[s]);
		efficiencies.push_back(_hVetoT[s]);

		_hHitPattern[s]
				= new TH2F(("hHitPattern" + toString(sensorIds[s])).c_str(), "Integrated hit pattern;Pixel X;Pixel Y",168,0,168,168,0,168);
		globals.push_back(_hHitPattern[s]);

		_hThreshEfficacy[s]
				= new TH1F(("hThreshEfficacy"+ toString(sensorIds[s])).c_str(), "Threshold Efficacy;Hit count per pixel", 20000, 0, 20000);
		globals.push_back(_hThreshEfficacy[s]);

		_hHitsInBtByThreshold[s]
				= new TH2F(("hHitsInBtByThreshold" + toString(sensorIds[s])).c_str(), "Hits in BT by threshold;Threshold; NHits", _bins, _lowEdge, _highEdge, 100, 0, 100);
		globals.push_back(_hHitsInBtByThreshold[s]);

		_hTimeStampsBySensor[s] = new TH1F(("hTimeStampsBySensor"+ toString(_sensorIds[s])).c_str(), ("hTimeStampsBySensor"+ toString(s)).c_str(), 40, 0, 40);
		globals.push_back(_hTimeStampsBySensor[s]);

		_hEffAndThreshSamp[s] = new TProfile(("hEffAndThreshSamp" + toString(sensorIds[s])).c_str(), ("hEffAndThreshSamp" + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		_hEffAndThreshShap[s] = new TProfile(("hEffAndThreshShap" + toString(sensorIds[s])).c_str(), ("hEffAndThreshShap" + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		_hEffAndThreshOverall[s] = new TProfile(("pEffAndThreshOverall" + toString(sensorIds[s])).c_str(), ("Overall efficiency, sensor " + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		efficiencies.push_back(_hEffAndThreshShap[s]);
		efficiencies.push_back(_hEffAndThreshSamp[s]);
		efficiencies.push_back(_hEffAndThreshOverall[s]);

		_hEffAndThreshSampDispT[s] = new TProfile(("hEffAndThreshSampDispT" + toString(sensorIds[s])).c_str(), ("hEffAndThreshSampDispT" + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		_hEffAndThreshShapDispT[s] = new TProfile(("hEffAndThreshShapDispT" + toString(sensorIds[s])).c_str(), ("hEffAndThreshShapDispT" + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		_hEffAndThreshOverallDispT[s] = new TProfile(("pEffAndThreshOverallDispT" + toString(sensorIds[s])).c_str(), ("Overall efficiency, sensor " + toString(sensorIds[s])).c_str(), _bins, _lowEdge, _highEdge);
		efficienciesDispT.push_back(_hEffAndThreshShapDispT[s]);
		efficienciesDispT.push_back(_hEffAndThreshSampDispT[s]);
		efficienciesDispT.push_back(_hEffAndThreshOverallDispT[s]);

		_hNHitsInBunchtrain[s] = new TH1F(("hNHitsInBunchtrain" + toString(sensorIds[s])).c_str(), ("hNHitsInBunchTrain" + toString(sensorIds[s])).c_str(), 1000000, 0, 1000000);
		dqm.push_back(_hNHitsInBunchtrain[s]);

		for (unsigned r(0); r < 4; r++) {

			_hHits[s][r] = new TH2F(("_hHitsRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), ("hHitsRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), _bins, _lowEdge, _highEdge, 400, 0, 400);
			hitsByRegion.push_back(_hHits[s][r]);

			_hWordHits[s][r] = new TH2F(("_hWordHitsRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), ("hWordHitsRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), _bins, _lowEdge, _highEdge, 400, 0, 4000);
			wordsByRegion.push_back(_hWordHits[s][r]);

			_hGroupsAndRows[s][r] = new TH2F(("_hGroupsAndRowsS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), ("hGroupsAndRowsS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), 7, -0.5, 6.5, 168, 0, 168);
			groupsAndRows.push_back(_hGroupsAndRows[s][r]);

			_hHitsByRegion[s][r] = new TProfile(("_hHitsByRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), ("_hHitsByRegionS" + toString(sensorIds[s]) + "R" + toString(r)).c_str(), _bins, _lowEdge, _highEdge);
			hitsByRegion.push_back(_hHitsByRegion[s][r]);

			_hTimestampsByThreshold[s][r] = new TH2F(("hTimestampsByThreshS"+ toString(sensorIds[s]) + "R" + toString(r)).c_str(), ("hTimestampsByThreshS"+ toString(sensorIds[s]) + "R" + toString(r)).c_str(), _bins, _lowEdge, _highEdge, 200, 0, 200);
			hitsByRegion.push_back(_hTimestampsByThreshold[s][r]);
		}
	}

	_hPmtTags1 = new TH1F("hPmtTags1", "hPmtTags1;Timestamp", 8192, 0, 8192);
	_hPmtTags2 = new TH1F("hPmtTags2", "hPmtTags2;Timestamp", 8192, 0, 8192);
	_hPmtCoincidence = new TH1F("hPmtCoincidence", "hPmtCoincidence;Timestamp", 8192, 0, 8192);
	pmt.push_back(_hPmtTags1);
	pmt.push_back(_hPmtTags2);
	pmt.push_back(_hPmtCoincidence);

	_hNHitsAndDispersion
			= new TH2F("hNHitsAndDispersion", "hNHitsAndDispersion;Dispersion;N_{hits}", 100, 0, 20, 15, 1, 16);
	_hNClustersInT = new TH1F("hNClustersInT", "hNClustersInT;Timestamp", 1024, 0, 8192);
	clustering.push_back(_hNHitsAndDispersion);
	clustering.push_back(_hNClustersInT);

	std::cout << "\tFinishing up in " << __PRETTY_FUNCTION__ << std::endl;
}

void MpsAnalysis::writeRootFile() {
	std::cout << __PRETTY_FUNCTION__ << std::endl;
	_rootFile->cd("/");
	for (tObjs::iterator it = globals.begin(); it != globals.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("hitsByRegion")->cd();
	for (tObjs::iterator it = hitsByRegion.begin(); it != hitsByRegion.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("wordsByRegion")->cd();
	for (tObjs::iterator it = wordsByRegion.begin(); it != wordsByRegion.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("groupsAndRows")->cd();
	for (tObjs::iterator it = groupsAndRows.begin(); it != groupsAndRows.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("efficiencies")->cd();
	for (tObjs::iterator it = efficiencies.begin(); it != efficiencies.end(); it++) {
		(*it)->Write();
	}

	_rootFile->mkdir("dispTeff")->cd();
	for (tObjs::iterator it = efficienciesDispT.begin(); it
			!= efficienciesDispT.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("dqm")->cd();
	for (tObjs::iterator it = dqm.begin(); it!= dqm.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("pmt")->cd();
	for (tObjs::iterator it = pmt.begin(); it!= pmt.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("clustering")->cd();
	for (tObjs::iterator it = clustering.begin(); it!= clustering.end(); it++) {
		(*it)->Write();
	}

	_rootFile->cd("/");
	_rootFile->mkdir("tracking")->cd();
	for (tObjs::iterator it = tracking.begin(); it!= tracking.end(); it++) {
		(*it)->Write();
	}
}