// -------------------------------------------------------------------------
// -----                 PndDpmDirect header file                  -----
// -----               Created 11/04/08  by M.Al-Turany              -----
// -------------------------------------------------------------------------

/** PndDpmDirect.h
 *@author M.Al-Turany <m.al-turany@gsi.de>
 *
 The PndDpmDirect generates DPM event using the DPM fortran code
 and inserts the tracks into the PndStack via the FairPrimaryGenerator.
 Derived from FairGenerator.
 *
 Filter for particle multiplicities added
 Minimum and Maximum multiplicities for each group id
 multiple conditions for the same PDG code added
 Martin dot J dot Galuska at physik dot uni minus giessen dot de
 **/

#ifndef PND_DPMDIRECT_H
#define PND_DPMDIRECT_H

#include "FairGenerator.h"
#include "TF1.h"
#include<map>
#include<vector>

class TVirtualMCStack;
class FairPrimaryGenerator;


typedef std::map<int, std::vector<int> >  PndPdgGroupNr;
typedef std::map<int,std::vector<int> >::iterator PndPdgGroupNrIterator;

std::ostream& operator <<(std::ostream& os, const std::vector<int>& v);
std::ostream& operator <<(std::ostream& os, const std::map<int, std::vector<int> >& pdgGroupNr);



class PndDpmDirect : public FairGenerator
{
public:

	/** Default constructor (should not be used) **/
	PndDpmDirect();


	/** Standard constructor
	 *  @param Mom in GeV/C
	 *  @param Mode = 0. - No elastic scattering, only inelastic
	 *  @param Mode = 1. - Elastic and inelastic interactions
	 *  @param Mode = 2. - Only elastic scattering, no inelastic one
	 **/
	PndDpmDirect(Double_t Mom, Int_t Mode, Long_t Seed = -1, Double_t ThtMin = -100.);
	PndDpmDirect(Double_t Mom, Int_t Mode, Double_t Rsigma, TF1* DensityFunction, Long_t Seed = -1, Double_t ThtMin=0.001);

	/** Destructor **/
	virtual ~PndDpmDirect();

	/** Generate one event using DPM
	 ** @param primGen  pointer to the FairPrimaryGenerator
	 **/
	virtual Bool_t ReadEvent(FairPrimaryGenerator* primGen);


	Bool_t AddFilterMinMax( int min, int max, int pdgCode1, int pdgCode2=0, int pdgCode3=0, int pdgCode4=0, int pdgCode5=0, int pdgCode6=0, int pdgCode7=0, int pdgCode8=0 )
	{
		// Use for grouping up to 8 pdgCodes into 1 groupId
		// min defines how many particles you want in your events AT LEAST
		// max defines how many particles you want in your events AT MOST

		// the min and max numbers are used for all particles with pdgCodes in the vector pdgCodes
		// this number is used for all particles with pdgCodes in the variables pdgCode1 to pdgCode8
		// pdgCodes groups different pdgCodes as one groupId (basically pretending that the listed pdgCodes are indistinguishable)


		// build a vector and call the general AddFilter
		std::vector<int> pdgCodes;
		pdgCodes.push_back(pdgCode1);
		pdgCodes.push_back(pdgCode2);
		pdgCodes.push_back(pdgCode3);
		pdgCodes.push_back(pdgCode4);
		pdgCodes.push_back(pdgCode5);
		pdgCodes.push_back(pdgCode6);
		pdgCodes.push_back(pdgCode7);
		pdgCodes.push_back(pdgCode8);

		return AddFilterMinMax( min, max, pdgCodes );
	}


	Bool_t AddFilterMinMax( int min, int max, std::vector<int> &pdgCodes )
	{
		// use if you want to group more than 8 pdgCodes into 1 groupId
		// min defines how many particles you want in your events AT LEAST
		// max defines how many particles you want in your events AT MOST

		// the min and max numbers are used for all particles with pdgCodes in the vector pdgCodes
		// pdgCodes groups different pdgCodes as one groupId (basically pretending that the listed pdgCodes are indistinguishable)
		// returns kTRUE if the filter could at least partially be added, otherwise returns kFALSE

		// 0 entries in pdgCodes will be skipped when setting the constraints


		// was adding of the filter successfully
		Bool_t filterAdded = kFALSE;

		if (fVerbose > 0){
			std::cout << "fPdgGroupNr: " << fPdgGroupNr << "\n";
			std::cout << "fPartCountsMin: " << fPartCountsMin << "\n";
			std::cout << "fPartCountsMax: " << fPartCountsMax << "\n";
			std::cout << "\n";
		}


		// check if filter can be applied (make sure that min and max make sense)
		if ( min < 1 ){
			std::cout << std::cout << "\n\n\n  -WARNING from PndDpmDirect: Filter could not be added. You are trying to request that your events should have at least a number <= 0 of some particles. That makes no sense. Check your AddFilterMinMax calls!\n\n\n";
			return kFALSE;
		}
		if ( max < min ){
			std::cout << std::cout << "\n\n\n  -WARNING from PndDpmDirect: Filter could not be added. You are trying to request that your events should have at most a number of some particles which is less than the minimum number that you request. That makes no sense. Check your AddFilterMinMax calls!\n\n\n";
			return kFALSE;
		}




		// find the groupId for the pdgCodes
		int groupId = fPartCountsMin.size();

		// check consistency of minimum and maximum multiplicity vector
		if ( fPartCountsMin.size() != fPartCountsMax.size() ){
			std::cout << std::cout << "\n\n\n  -FATAL ERROR from PndDpmDirect: fPartCountsMin and fPartCountsMax have different sizes. That should not happen! Filter is not added!\n\n\n";
			return kFALSE;
		}


		// for every entry in pdgCodes
		// -- save a new entry in the map fPdgGroupNr if the entry is not already in map
		// -- add a new groupId to the vector in fPdgGroupNr if the entry already exists in the map
		typedef std::pair<int,std::vector<int> > PndPdgGroupNrPair;

		for (int iPdgCodes = 0; iPdgCodes < pdgCodes.size(); ++iPdgCodes){
			// skip 0 entries in pdgCodes
			if ( 0 != pdgCodes[iPdgCodes] ) {
				// check if map has already an entry for the pdgCode
				PndPdgGroupNrIterator finder = fPdgGroupNr.find(pdgCodes[iPdgCodes]);
				if ( fPdgGroupNr.end() != finder){
					// pdgCode is already in map, add groupId to vector

					std::vector<int> & groupIdVector = fPdgGroupNr[pdgCodes[iPdgCodes]];
					groupIdVector.push_back(groupId);
				}
				else {
					// pdgCode is not in map yet, create new entry of int and vector<int>

					std::vector<int> groupIdVector;
					groupIdVector.push_back(groupId);
					fPdgGroupNr.insert(PndPdgGroupNrPair(pdgCodes[iPdgCodes],groupIdVector));
				}

				filterAdded = kTRUE;
			}
		}

		if ( kTRUE == filterAdded ){
			// set the minimum and maximum number of particles for the current groupId
			fPartCountsMin.push_back(min);
			fPartCountsMax.push_back(max);
		} else {
			std::cout << "\n\n\n  -WARNING from PndDpmDirect: Tried to add filter, but that was not possible as all pdg codes were 0!\n\n\n";
			return kFALSE;
		}

		// turn filter on
		fFilterEvents = kTRUE;




		if (fVerbose > 0){
			std::cout << "PndDpmDirect: After adding Filter:\n";
			std::cout << "fPdgGroupNr: " << fPdgGroupNr << "\n";
			std::cout << "fPartCountsMin: " << fPartCountsMin << "\n";
			std::cout << "fPartCountsMax: " << fPartCountsMax << "\n";
			std::cout << "\n";
		}

		return kTRUE;
	}

	void SetFilterMaxTries(int maxTries=99999)
	{
		if ( maxTries > 0 ){
			fFilterMaxTries = maxTries;
			std::cout << "PndDpmDirect: maxTries is now set to " << fFilterMaxTries << "\n";
		} else {
			std::cout << "\n\n\n  -WARNING from PndDpmDirect: maxTries must be a positive number! Check your SetFilterMaxTries call!\n\n\n";
		}
	}

	int GetNumberOfSimulatedEvents(){
		return fSimulatedEvents;
	}

	int GetNumberOfFilterFailedEvents(){
		// this should return 0 if everything works fine
		// if it returns a value >0 it means that fFilterMaxTries is set too low
		// if it returns a value which is equal to the number of events that you requested, it means that either fFilterMaxTries is set way too low
		// or that the generator does not create such events that you are interested in
		// or that your conditions cannot be satisfied at all (logical error)
		return fFailedEvents;
	}


	void SetVerbose(int verbose=11)
	{
		fVerbose = verbose;
	}


private:


	Bool_t fFilterEvents; // Filter the events produced by dpm iif kTRUE
	PndPdgGroupNr fPdgGroupNr; // first input = PID (as used by dpm generator) second input = "vector of group number" // vector is not implemented yet
	// particles with different pid can be counted as the same by assigning them the same group number
	// group numbers have to start with 0 and be consecutive as they are used as the index for fPartCounts
	std::vector<int> fPartCountsMin; // vector containing the minimum number of particles (according to "group number") for acceptable events
	std::vector<int> fPartCountsMax; // vector containing the maximum number of particles (according to "group number") for acceptable events
	int fSimulatedEvents;
	int fFilterMaxTries;
	int fFailedEvents;

	Int_t fVerbose;

	/**
	 * P_lab(GeV/c)
	 */
	float fMom;
	/**
	 * 0. - No elastic scattering, only inelastic
	 * 1. - Elastic and inelastic interactions
	 * 2. - Only elastic scattering, no inelastic one
	 */
	float  fMode;
	double fSeed;
	int    fGasmode;
	double fRsigma;
	float fThtMin;

	TF1*   fDensityFunction; //!

	ClassDef(PndDpmDirect,1);

};

#endif