// -------------------------------------------------------------------------
// -----                PndDpmDirect source file                  -----
// -----                                                          -----
// -------------------------------------------------------------------------


#include <iostream>
#include "TClonesArray.h"
#include "TFile.h"
#include "TLorentzVector.h"
#include "TTree.h"
#include "TVector3.h"
#include "TParticle.h"
#include "PndDpmDirect.h"
#include "FairPrimaryGenerator.h"
#include "TRandom.h"

using namespace std;

#include "cfortran.h"
extern "C" {
// COMMON /LUJETS/ N,K(1000,2),P(1000,5)
// n   - number of produced particles,
// k[] - Pythia particle identifiers
// p[] - kinematical characteristics of particles


typedef struct {
	int n, k[2000];
	float p[5000];
} lujets;

#define LUJETS COMMON_BLOCK(LUJETS,lujets)
COMMON_BLOCK_DEF(lujets,LUJETS);

}


std::ostream& operator <<(std::ostream& os, const std::vector<int>& v)
{
	for (int iVec = 0; iVec < v.size(); ++iVec){
		os << v[iVec] << ", ";
	}
	os << "])  ";
	return os;
}


std::ostream& operator <<(std::ostream& os, const std::map<int, std::vector<int> >& pdgGroupNr)
{
	for (std::map<int,std::vector<int> >::const_iterator iPdgGroupNr = pdgGroupNr.begin(); iPdgGroupNr != pdgGroupNr.end(); ++iPdgGroupNr){
		os << "("<< iPdgGroupNr->first << ", [";
		os << iPdgGroupNr->second;
		os << "])  ";
	}
	return os;
}


extern "C" int init1_(float* Plab, double* seed, float* Elastic, float* tetmin );//install DPM
extern "C" int dpm_gen_(float* Generator, double* seed ); //to generate events

TF1 	 * fDensityFunction;


// -----   Default constructor   ------------------------------------------
PndDpmDirect::PndDpmDirect() : fFilterEvents(kFALSE), fSimulatedEvents(0), fVerbose(0), fFailedEvents(0), fFilterMaxTries(9999) {

}
// ------------------------------------------------------------------------

// -----   Standard constructor   -----------------------------------------

PndDpmDirect::PndDpmDirect(Double_t Mom, Int_t Mode, Long_t Seed, Double_t ThtMin) : fFilterEvents(kFALSE), fSimulatedEvents(0), fVerbose(0), fFailedEvents(0), fFilterMaxTries(9999)  {
	fMom    = Mom;
	fMode   = Mode;
	fSeed   = Seed;
	fThtMin = ThtMin;


	if ( ThtMin < -90 ){
		// Calculate ThtMin first. For this we make a cut-off on the value of -t of 1e-2 GeV^2 (~100 MeV/c momentum)
		// This estimated from a parametrization found in thesis of Thomas Wuerschig (figure 6.4, page 121):
		// Roughly: 0.4 deg at 15 GeV/c and 4 deg at 1.5 GeV/c, linear interpolation in double log-scale.
		//
		Double_t logangle = TMath::Log(0.4)+(TMath::Log(15.)-TMath::Log(Mom))*(TMath::Log(4)-TMath::Log(0.4))/(TMath::Log(15)-TMath::Log(1.5));
		ThtMin = TMath::Exp(logangle);
	}

	if (fSeed < 0)
	{
		Long_t iSeed = gRandom->GetSeed();
		int a = iSeed/100000;
		fSeed=iSeed - a*100000 + a/100000.;
	}

	fGasmode = 0;
	fRsigma = 0.;
	fDensityFunction = new TF1();

	cout << "<I> PndDpmDirect initialization" << endl;
	cout << "<I> Momentum  = " << fMom << endl;
	cout << "<I> Seed      = " << fSeed << endl;
	cout << "<I> Mode      = " << fMode << endl;
	cout << "<I> Theta min = " << fThtMin<<endl;

	init1_(&fMom,&fSeed,&fMode, &fThtMin);     // init the DPM generator
}
// ------------------------------------------------------------------------

// -----   Gas mode constructor   -----------------------------------------
PndDpmDirect::PndDpmDirect(Double_t Mom, Int_t Mode, Double_t Rsigma, TF1* DensityFunction, Long_t Seed, Double_t ThtMin) : fFilterEvents(kFALSE), fSimulatedEvents(0), fVerbose(0), fFailedEvents(0), fFilterMaxTries(9999)  {
	fMom    = Mom;
	fMode   = Mode;
	fSeed   = Seed;
	fThtMin = ThtMin;

	if (fSeed < 0)
	{
		Long_t iSeed = gRandom->GetSeed();
		int a = iSeed/100000;
		fSeed=iSeed - a*100000 + a/100000.;
	}

	fGasmode = 1;
	fRsigma = Rsigma;
	fDensityFunction = DensityFunction;

	cout << "<I> PndDpmDirect initialization" << endl;
	cout << "<I> Momentum  = " << fMom << endl;
	cout << "<I> Seed      = " << fSeed << endl;
	cout << "<I> Mode      = " << fMode << endl;
	cout << "<I> Gasmode   = " << fGasmode << endl;
	cout << "<I> Theta min = " << fThtMin<<endl;

	init1_(&fMom,&fSeed,&fMode, &fThtMin);     // init the DPM generator
}
// ------------------------------------------------------------------------

// -----   Destructor   ---------------------------------------------------
PndDpmDirect::~PndDpmDirect() {

}
// ------------------------------------------------------------------------

// -----   Public method ReadEvent   --------------------------------------
Bool_t PndDpmDirect::ReadEvent(FairPrimaryGenerator* primGen) {
	int  npart, i;
	Double_t fX, fY, fZ, radius;
	double Px[1000],Py[1000],Pz[1000],E[1000],Pm[1000],Wh[1000];
	int Id[1000];

	float Generator=0.; //  Format in which events are produced (0=pythia, 1=pluto)

	Double_t weight = 1.0;
	Int_t activeCnt=0;



	if (fVerbose > 3){
		std::cout << "PndDpmDirect: Beginning of ReadEvent\n";
		std::cout << "fPdgGroupNr: " << fPdgGroupNr << "\n";
		std::cout << "fPartCountsMin: " << fPartCountsMin << "\n";
		std::cout << "fPartCountsMax: " << fPartCountsMax << "\n";
		std::cout << "\n";
	}




	{
		// check if event is acceptable
		Bool_t evtOk;

		// count how often dpm is run to find a suitable event
		Int_t iTry = 0;

		do{
			// run generator
			dpm_gen_(&Generator, &fSeed);

			// get all produced particles
			npart = lujets_.n;

			// check the produced particles if you want to filter the generated events
			if ( kTRUE == fFilterEvents ) {
				// initialise a 0 counter of same size as fPartCountsMin / Max
				std::vector<int> counter;
				for (int iCounter = 0; iCounter < fPartCountsMin.size(); ++iCounter){
					counter.push_back(0);
				}
				// counter counts up, event is accepted if all entries in counter are between the according entries in fPartCountsMin and fPartCountsMax
				if (fVerbose >11){
					cout << "counter after initialization: " << counter << "\n";
				}


				// loop the produced particles and count up the needed "groupIDs" using the PIDs specified in fPdgGroupNr
				for (Int_t iPart= 0; iPart< npart; ++iPart) {

					// get pdg code for particle iPart
					int partId = lujets_.k[iPart+1000];

					// check if the pdg code is contained in some groupId
					PndPdgGroupNrIterator it = fPdgGroupNr.find(partId);
					if ( it != fPdgGroupNr.end() ){
						// TODO: Check momentum and angles... for all particles
						// count up all groupIds which satisfy the momentum/angle conditions (in a loop over the vector entries)
						std::vector<int> groupIdVector = it->second;
						for (int iGroupId = 0; iGroupId < groupIdVector.size(); ++iGroupId){
							++counter[groupIdVector[iGroupId]];
						}
					}

					if (fVerbose >11){
						cout << "iTry == " << iTry << "\n";
						cout << "counter: " << counter << "\n";
					}

				}


				// reset evtOk to kTRUE for each new event
				evtOk = kTRUE;

				if (fVerbose >9){
					cout << "Check if event is acceptable \n";
					cout << "evtOk is " << evtOk << "\n";
					cout << "counter: " << counter << "\n";
					cout <<  "fPartCountsMin: " << fPartCountsMin << "\n";
					cout <<  "fPartCountsMax: " << fPartCountsMax << "\n";
				}

				// check if event is acceptable
				// event is acceptable if
				// counter[iCounter] >= fPartCountsMin[iCounter]
				// and
				// counter[iCounter] <= fPartCountsMax[iCounter]
				for (int iCounter = 0; iCounter < counter.size(); ++iCounter){
					if ( counter[iCounter] < fPartCountsMin[iCounter] ){
						evtOk = kFALSE;
						if (fVerbose >9){
							cout << "Set evtOk to " << evtOk << " because of counter[iCounter] < fPartCountsMin[iCounter] for iCounter == " << iCounter << "\n";
						}
					}
					if ( counter[iCounter] > fPartCountsMax[iCounter] ){
						evtOk = kFALSE;
						if (fVerbose >9){
							cout << "Set evtOk to " << evtOk << " because of counter[iCounter] > fPartCountsMax[iCounter] for iCounter == " << iCounter << "\n";
						}
					}
				}
			}

			if ( kTRUE == evtOk ) {
				if (fVerbose >5){
					cout << "Event is ok. Needed " << iTry + 1 << " tries to succeed!\n";
				}

			}
		} while ( ++iTry<fFilterMaxTries && !evtOk );
		fSimulatedEvents+=iTry;
		if ( !evtOk ){
			++fFailedEvents;
			cout << "\n -E PndDpmDirect: No event was found within " << fFilterMaxTries << " which satisfies your event filter. Try increasing the max. number of tries or change your filter\n\n";
		}
		if (fVerbose > 4){
			cout << fSimulatedEvents << " events simulated until I found a good one or gave up.\n";
			cout << fFailedEvents << " unsuccessful attempts to find an event that suits your filters\n";
		}
	}


	for (i= 0; i< npart; ++i) {

		Id[i]=lujets_.k[i+1000];
		Px[i]=lujets_.p[i];
		Py[i]=lujets_.p[i+1000];
		Pz[i]=lujets_.p[i+2000];
		Pm[i]=lujets_.p[i+4000];
		E[i]=lujets_.p[i+3000];
		Wh[i]=1.0;

		/* Check if fGasmode is set */
		fX = 0.;
		fY = 0.;
		fZ = 0.;
		if (fGasmode == 1) {

			// define position of track start
			// Random 2D point in a circle of radius r (simple beamprofile)
			radius = gRandom->Gaus(0,fRsigma);
			gRandom->Circle(fX, fY, radius);

			// calculate fZ according to some (probability) density function of the gas
			fZ=fDensityFunction->GetRandom();

		}

		// add track
		//printf("- I -: new particle at: %f, %f, %f ...\n", fX, fY, fZ);
		primGen->AddTrack(Id[i], Px[i], Py[i], Pz[i], fX, fY, fZ);

	}

	return kTRUE;

}
// ------------------------------------------------------------------------

ClassImp(PndDpmDirect)