//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
//olgusha.net// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//---------------------------------------------------------------------------
//
// ClassName:   Pixel - Charge transport and sharing
//
// Author:      M. Esposito 20.03.13
//
//----------------------------------------------------------------------------
//
#include "Pixel.hh"
#include "Transport.hh"
///#include "Histo.hh"
#include "Randomize.hh"
#include "Pixel.hh"
#include "math.h"
#include <time.h>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/storage.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/lognormal_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <sstream>

#include <map>

#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"

using namespace boost::numeric::ublas;

Pixel::Pixel()
{
	///histo = new Histo();
	histoManager = HistoManager::GetHistoManager();
	
	pixelMessenger = new PixelMessenger(this);

	// default values
	length = 100*mm;
	pitch = 200*um;
	threshold=500;
}

Pixel::~Pixel(){	
	///delete histo;
	if(histoManager){
		delete histoManager;
	}
}

void Pixel::Sample(std::vector <G4double> xHit,std::vector <G4double> yHit, G4int motherCopyno, G4String RTtype,G4int number, G4int PDGid)
{	

///	G4cout << "Pixel::Sample::Test" << G4endl;

// *************************************************************************************************************
//********************* Detector parameters ********************************************************************
//**************************************************************************************************************

	if (RTtype=="CMOS"){
		pitch=pitch;/// 0.05*mm;
	}
	else{
		pitch= 0.100*mm;
	}
	
	///G4cout << "length = " << length << ", pitch = " << pitch << ", threshold = " << threshold << G4endl;

	unsigned noPixel=length/pitch;
	
	// get the number of electrons as passed by transport
	G4double electrons=xHit.size();

	//Parameters for DynAMITE: not currently used
	G4double gainM=59; //e/DN mean of the lognormal distribution
	G4double gainS=21; //e/DN std of the lognormal distribution
	//G4double noiseM=262; //e mean of the lognormal distribution
	G4double noiseM=698; //e mean of the lognormal distribution 8 LSB noise
	G4double noiseS=168; //e std of the lognormal distribution
	
//**************************************************************************************************************
//**************************************************************************************************************

// *************************************************************************************************************
//********************* Matrices *******************************************************************************
//**************************************************************************************************************
/*	//Hit: Nx2 matrix with pixel coordinates
	matrix<double> Hit(electrons,2);
	// Frame: 2Nx2N oversampled matrix counting hit per pixel
	matrix<double> Frame(noPixel+1,noPixel+1);
	// SharedFrame: 2Nx2N oversampled matrix with charged shared in the neighbour pixels
	matrix<double> SharedFrame(noPixel*2,noPixel*2);
	//ResFrame: NxN matrix resampled after charge sharing
	matrix<double> ResFrame(noPixel,noPixel);
	//DNFrame: NxN matrix resampled after charge sharing and convert to DN
	matrix<double> DNFrame(noPixel,noPixel);
//**************************************************************************************************************
//**************************************************************************************************************

// *************************************************************************************************************
//********************* Random Engine **************************************************************************
//**************************************************************************************************************
	typedef boost::mt19937 ENG;    // Mersenne Twister
    typedef boost::lognormal_distribution<double> DISTg;   // LogNormal Distribution Gain
    typedef boost::variate_generator<ENG,DISTg> GENg;    // Variate generator Gain
 	typedef boost::lognormal_distribution<double> DISTn;   // LogNormal Distribution Noise
    typedef boost::variate_generator<ENG,DISTn> GENn;    // Variate generator Noise
    ENG  eng;
    eng.seed(time(0));
    DISTg distG(gainM,gainS);
    GENg  genG(eng,distG);
    DISTn distN(noiseM,noiseS);
    GENn  genN(eng,distN);
 //**************************************************************************************************************
//**************************************************************************************************************
    G4int pixel=0;
    G4int x=0;
    G4int y=0;

		// set the pixel values all to zero
    for (unsigned i=0; i<noPixel; i ++){
    	for(unsigned j=0; j<noPixel; j++){
			Frame(i,j)=0;
		}
    }

	
	// set the xy position of all electrons to be zero i.e reset matrix
	// this is not actually needed as the next loop sets positions anyway
    for(unsigned i=0; i< electrons; i++){
		Hit(i,0)=0;
		Hit(i,1)=0;				
	}

	//Filling Hit with Pixel coordinates of all the electrons
	for(unsigned i=0; i<electrons; i++){

		Hit(i,0)=floor((xHit[i]+(length/2))/(pitch));
		Hit(i,1)=floor((yHit[i]+(length/2))/(pitch));
	}
	
	//if we have some electrons loop through and count up how many in each pixel
	if(electrons){
		// //Filling Frame
		// loop throught the hit matrix
		for (unsigned i=0; i< Hit.size1();i++){
			// if the hit is within the sensor then use
			if(Hit(i,0)<noPixel && Hit(i,1)<noPixel && Hit(i,0)>=0 && Hit(i,1)>=0  ){
				//incriment the number of electrons in the location by 1 from the position of the hit
				Frame(Hit(i,0),Hit(i,1))=Frame(Hit(i,0),Hit(i,1))+1;
			}
		}

	
	
	//passing output to histo for write on file
	//loop through all of this pixels and see which ones contain enough electrons above threshold
	// i think this is bottle neck as looping 2000x2000 each time for sparse data
	for (unsigned i=0; i<noPixel; i ++){
    	for(unsigned j=0; j<noPixel; j++){
			if(Frame(i,j)>threshold){
///				G4cout << Frame(i,j) << G4endl;
// 				data.push_back(number);
// 				data.push_back(PDGid);
// 				data.push_back(motherCopyno);
// 				data.push_back(i);
// 				data.push_back(j);
// 				data.push_back(Frame(i,j));
// 				G4cout << Frame(i,j) << G4endl;
// 	///			histo->FillBuffer(data);
// 				data.clear();
				histoManager->AddRTEvent(number, PDGid, motherCopyno, i, j, Frame(i,j));

			}	
		}
   }
}  */ 
// // // 
// // //     //This is to implement charge division between adjent pixel and conversion in DN
// // //    //Not recommended for Priapus yet as performance not known
// // // 
// // // 	// //Filling SharedFrame
// // // 	// for(unsigned i=1; i<Frame.size1(); i++){
// // // 	// 	for (unsigned j=1;j<Frame.size2();j++){
// // // 	// 		if(Frame(i,j)){
// // // 	// 			SharedFrame(i,j)=SharedFrame(i,j)+(Frame(i,j)/4);
// // // 	// 			SharedFrame(i-1,j)=SharedFrame(i-1,j)+(Frame(i,j)/4);
// // // 	// 			SharedFrame(i-1,j-1)=SharedFrame(i-1,j-1)+(Frame(i,j)/4);
// // // 	// 			SharedFrame(i,j+1)=SharedFrame(i,j+1)+(Frame(i,j)/4);
// // // 	// 		}		
// // // 	// 	}
// // // 	// }
// // // 	// G4double sumSh=0;
// // // 	// G4double sumRes=0;
// // // 	// G4double sumDN=0;
// // // 	// G4double sumDNnoise=0;
// // // 	// G4int count=0;
// // // 	// G4double gainCurr=0;
// // // 	// G4double noiseCurr=0;
// // // 	// // Resampling shared Frame
// // // 	// unsigned m=0;
// // // 	// unsigned n=0;
// // // 	// for(unsigned i=1; i<noPixel; i++){
// // // 	// 	m=i*2;
// // // 	// 	for (unsigned j=1;j<noPixel;j++){
// // // 	// 		n=j*2;
// // // 	// 		ResFrame(i,j)=Frame(m,n)+Frame(m,n+1)+Frame(m+1,n)+Frame(m+1,n+1);
// // // 	// 		sumRes=sumRes+ResFrame(i,j);
// // // 	// 		if(ResFrame(i,j)){
// // // 	// 			//G4cout << i << " " << j << G4endl;
// // // 	// 			gainCurr=genG();
// // // 	// 			DNFrame(i,j)=ResFrame(i,j)/gainCurr;
// // // 	// 			count=count+1;
// // // 	// 			sumDN=sumDN+DNFrame(i,j);
// // // 	// 			noiseCurr=genN();
// // // 	// 			//G4cout << "Noise " << noiseCurr << " Gain " << gainCurr << G4endl;
// // // 	// 			DNFrame(i,j)=(ResFrame(i,j)-noiseCurr)/gainCurr;
// // // 	// 			sumDNnoise=sumDNnoise+DNFrame(i,j);
// // // 	// 		}
// // // 	// 	}
// // // 	// }
// // // 	//G4cout << sumDN << G4endl;
// // // 
// // // 	
// // // 	
// // // 
// // // 	}
// // // 	

	// this condenses the above 3 for loops into ones
	// speeds this section of code up x3

	if(electrons)
	{
		std::map<std::pair<unsigned,unsigned>, unsigned> Frame;
		for(unsigned i=0; i<electrons; i++)
		{
			unsigned xpix = floor((xHit[i]+(length/2))/(pitch));
			unsigned ypix = floor((yHit[i]+(length/2))/(pitch));
			
			std::pair<unsigned,unsigned> hit(xpix,ypix);
			
			if(xpix<noPixel && ypix<noPixel && xpix>=0 && ypix>=0)
			{
				//Frame(xpix,ypix)=Frame(xpix,ypix)+1; // eventually this will be changed to a map
				std::map<std::pair<unsigned,unsigned>, unsigned>::iterator it = Frame.find(hit);
				if(Frame.find(hit) != Frame.end()) // if already exists then add an electron to the frame(xpix,ypix)
					(*it).second++;
				else                               // if it does not exist add it to the map and set nelectrons = 1
					Frame[hit] = 1;
			}
		}
		
		// now loop through the map and check if above threshold
		std::map<std::pair<unsigned,unsigned>, unsigned>::iterator it;
		for(it=Frame.begin(); it!=Frame.end(); it++)
		{
			if( (*it).second > threshold )
			{
				if(histoManager){
				histoManager->AddRTEvent(number, PDGid, motherCopyno, (*it).first.first, (*it).first.second, (*it).second);
				}
				else{
					G4cout << number << " " <<  PDGid << " " << motherCopyno << " " << (*it).first.first << " " << (*it).first.second << " " << (*it).second << G4endl;
				}
			}
		}
	}

}