#include "TString.h"
#include "TObjArray.h"
#include "TObjString.h"


#include "TFile.h"
#include "TNtuple.h"
#include "TSystem.h"
#include "TRandom.h"

#include "hphysicsconstants.h"
#include "hseed.h"

#include "PParticle.h"
#include "PChannel.h"
#include "PFixedDecay.h"
#include "PFireball.h"
#include "PDoubleFireball.h"
#include "PReaction.h"
#include "PEmbeddedParticles.h"
#include "PHGeantOutput.h"
#include "PVertexFile.h"

#include <fstream>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>


//----------------------------------------------------------------
// DATABASE HELPERS
typedef struct {
  TString type;
  TString channel;
  Float_t T1;
  Float_t T2;
  Float_t fracT1;
  Float_t blast;
  Float_t A2;
  Float_t A4;
  Float_t V1;
  Float_t V2;
  Float_t eta;

  PFireball* source;
  vector<PChannel*>  vch;         // list of channels
  vector<PParticle*> vm;          // mother of decays
  vector<vector<PParticle*> > vd; // daughters of decays
  PChannel** cfinal;
  void clear(){
      type = "";
      channel = "";
      T1   = 0;
      T2   = 0;
      fracT1 = 0;
      blast= 0;
      A2   = 0;
      A4   = 0;
      V1   = 0;
      V2   = 0;
      eta  = 0;
      source = 0;

      cfinal = 0;
      vch.clear();
      vm.clear();
      for(UInt_t i = 0; i < vd.size(); i++) { vd[i].clear();}
  }

  void print(Int_t i=-1){

      if(i < 0) {
	  cout<<setw(20)<<type
	      <<" chan "  <<setw(30)<<channel
	      <<" T1: "   <<setw(10)<<T1
	      <<" T2: "   <<setw(10)<<T2
	      <<" fracT1: "<<setw(10)<<fracT1
	      <<" blast: "<<setw(10)<<blast
	      <<" A2: "   <<setw(10)<<A2
	      <<" A4: "   <<setw(10)<<A4
	      <<" V1: "   <<setw(10)<<V1
	      <<" V2: "   <<setw(10)<<V2
	      <<" eta: "  <<setw(10)<<eta
	      <<endl;
      } else {
	  cout<<setw(3)<<i<<" : "<<setw(20)<<type
	      <<" chan "  <<setw(30)<<channel
	      <<" T1: "   <<setw(10)<<T1
	      <<" T2: "   <<setw(10)<<T2
	      <<" fracT1: "<<setw(10)<<fracT1
	      <<" blast: "<<setw(10)<<blast
	      <<" A2: "   <<setw(10)<<A2
	      <<" A4: "   <<setw(10)<<A4
	      <<" V1: "   <<setw(10)<<V1
	      <<" V2: "   <<setw(10)<<V2
	      <<" eta: "  <<setw(10)<<eta
	      <<endl;
      }
  }

} thermaldecay;

struct database {
    TString beamtime;
    map<Float_t, vector<thermaldecay*> > mEtoDecay;

    void print(){
	cout<<"#######################################################################################"<<endl;
	cout<<"DATABASE Beamtime = "<<beamtime<<endl;

	map<Float_t,vector<thermaldecay*> >::iterator it;
        Int_t ct=1;
	for(it = mEtoDecay.begin(); it != mEtoDecay.end(); ++it ){
            cout<<"---------------------------------------------------------------------------------------"<<endl;
	    cout<<"SET "<<ct<<": Beamenergy = "<<it->first<<endl;
	    for(UInt_t i = 0 ; i < it->second.size(); i++){
               it->second[i]->print(1+i);
	    }
            ct++;
	}
	cout<<"---------------------------------------------------------------------------------------"<<endl;
	cout<<"#######################################################################################"<<endl;
    }

    thermaldecay* getDecayParams(Bool_t& ret,Float_t energy, TString type)
    {
	ret = kFALSE;
	map<Float_t,vector<thermaldecay*> >::iterator it;
	it = mEtoDecay.find(energy);
	if(it == mEtoDecay.end()) {
	    cout<<"ERROR : energy = "<<energy<<" not found in dta base"<<endl;
	    return 0;
	} else {
            for(UInt_t i = 0 ; i < it->second.size(); i++){
		if(it->second[i]->type == type){
                    ret = kTRUE;
		    return it->second[i];
		}
	    }
	    cout<<"ERROR : type = "<<type.Data()<<" not found in data base for energy = "<<energy<<endl;
            return 0;
	}
	return 0;
    }

    Int_t getParticle(TString list,vector<PParticle*>& vp,Bool_t print=kFALSE){

	vp.clear();
	TString tmp = list;
	TObjArray* a = tmp.Tokenize(",");
	if(a){
	    for(Int_t i = 0 ; i < a->GetEntries(); i++){
		TObjString* s = (TObjString*) a->At(i);
		TString t = s->GetString();
		vp.push_back(new PParticle(t.Data()));
		if(print) cout<<"getParticle() : "<<i<<" "<<t<<endl;
	    }

	    delete a;
	}
	return vp.size();
    }

    void parseChannel(Float_t Energy,thermaldecay* d)
    {

	TString channel = d->channel;
	channel.ReplaceAll(" ","");          // no spaces allowed in channel descriptor

	Ssiz_t i = channel.First(":");       // format examples:  "eta:e+,e-,..."  or  'pi0:dilepton[e+e-],g'
	TString mother = channel(0,i);       // main particle for the channel
	channel.Remove(0, i+1);


	//----------------------------------------------
	// looping the channel descriptor and defining the different decay levels
	TString l   [10]; // level of decay by appearance of '[' in channel (max 10), comma separated particles as string
	TString ml  [10]; // decay mother of level as string
	Int_t   ind [10];
	Int_t   indm[10];

	for(Int_t i = 0; i < 10; i++) { ind[i] = 0; indm[i] = 0; }

	vector<PParticle*> vl [10]; // PParticles of decay level
	Int_t level = 0;
	Int_t max   = 0;

	ml[0] = mother;                               // level 0 has the main particle as mother
	d->vm.push_back(new PParticle(ml[0].Data()));  // create a PParticle for the main particle

	for(Int_t i = 0; i < channel.Length(); i++){
	    Char_t c = channel[i];
	    if (c == '[' ) { // new decay
		level++;
		if(level > max ) max = level;

		Int_t j = i;  // loop back to last , or begin of string to find mother of decay
		for(; j > 0; j--){
		    Char_t ds = channel[j];
		    if (ds == ',' || ds == '[' ) {
			j++;
			break;
		    }
		}

		ml  [level] = channel(j,i-j);  // mother of this decay
		indm[level] = ind[level-1];    // remember the index of mother particle
		continue;
	    }

	    if (c == ']' ) { // end decay
		level--;
		continue;
	    }

	    if (c == ',') ind[level]++; // remember the index of mother particle
	    l[level] += channel[i];     // will end up as comma separated particle list per level
	}
	//----------------------------------------------


	//----------------------------------------------
	// setup Fireball or DoubleFireball source
	if(d->eta == 0){
	    d->source = new PFireball      ((Char_t*)ml[0].Data() ,Energy,d->T1,d->T2,d->fracT1,d->blast,       d->A2,d->A4,d->V1,d->V2);
	} else {
	    d->source = new PDoubleFireball((Char_t*)ml[0].Data() ,Energy,d->T1,d->T2,d->fracT1,d->blast,d->eta,d->A2,d->A4,d->V1,d->V2);
	}

	d->source->setTrueThermal(kTRUE);
	d->source->Print();

	PParticle *moth = d->vm[0];  // the mother of the decay is the main particle

	PParticle** smoth  = new PParticle* [2];  // connect mother to Fireball
	smoth[0] = d->source;
	smoth[1] = moth;

	PChannel* cmoth = new PChannel(smoth, 1, 1);  // create a channel for the Fireball
	d->vch.push_back(cmoth);                       // add to the list of channels
	//----------------------------------------------

	if( !(max == 0 && mother == l[0])) {  // not a fireball particle without decay inside PLUTO

	    //----------------------------------------------
	    // analyze the level strings to build the decay channels
	    // for PReaction
	    for(Int_t i = 0 ; i < max+1; i++){
		cout<<"decay level "<<i<<endl;
		getParticle(l[i],vl[i],kTRUE); // comma separated particle string -> PParticles
		d->vd.push_back(vl[i]);

		if(i > 0) { // find the mother of decay
		    d->vm.push_back(vl[i-1][indm[i]]) ;
		    TString n = HPhysicsConstants::pid(d->vm[i]->ID());
		    if( n != ml[i] ){
			cout<<"############################################################"<<endl;
			cout<<"ERROR : mismatch pid from index with stored name : pid from index = "<<HPhysicsConstants::pid(d->vm[i]->ID())<<" stored "<<ml[i]<<" index "<<indm[i]<<endl;
			cout<<"############################################################"<<endl;
		    }
		}

		//----------------------------------------------
		// build Channel for this decay level
		PParticle **s = new PParticle* [vl[i].size()+1];
		s[0] = d->vm[i] ; // mother of this decay
		for(UInt_t j = 0; j < vl[i].size(); j++){
		    s[j+1] = vl[i][j]; // add daughters
		}
		d->vch.push_back(new PChannel(s,vl[i].size(),1)); // add to channel list
		//----------------------------------------------
	    }
	}

	//----------------------------------------------
	// build the final channel from collect channels above
	cout<<"############################################################"<<endl;
	cout<<"CONFIG THERMAL PARTICLE : "<<d->type<<endl;
	cout<<"FINAL CHANNEL SETUP for : channel descriptor : "<<d->channel<<endl;
	d->cfinal = new PChannel* [d->vch.size()];
	for(UInt_t i = 0 ; i < d->vch.size(); i++){
	    d->cfinal[i] = d->vch[i];
	    d->cfinal[i] ->PrintReaction();
	}
	cout<<"############################################################"<<endl;
	//----------------------------------------------

    }

    Bool_t readDataBase(database& datBase, TString infile)
    {

	if(gSystem->AccessPathName(infile.Data()) != 0){

	    cout<<"############################################################"<<endl;
	    cout<<"readDataBase() : data base input file "<<infile<<" not found!"<<endl;
	    cout<<"############################################################"<<endl;
	    return kFALSE;
	}
	cout<<"############################################################"<<endl;
	cout<<"READING DATABASE INPUT FILE : "<<infile<<endl;
	cout<<"############################################################"<<endl;

	Char_t line[800];
	ifstream in;
	in.open(infile.Data(), std::ifstream::in);

	Float_t energy = 0;

	map<Float_t,vector<thermaldecay*> >::iterator it;

	while(in.good()){
	    in.getline(line,800);
	    TString tmp = line;
	    tmp.ReplaceAll(" ","");

	    if(tmp.Contains("//"))   continue;
	    if(tmp.Contains("#"))    continue;
	    if(tmp.CompareTo("")==0) continue;

	    if(tmp.Contains("beamtime")){
		tmp.ReplaceAll("beamtime=","");
		datBase.beamtime = tmp;
		continue;
	    }

	    if(tmp.Contains("beamenergy")){
		tmp.ReplaceAll("beamenergy=","");
		energy = tmp.Atof();
		continue;
	    }

	    thermaldecay* d = new thermaldecay;
	    d->clear();

	    TObjArray* a = tmp.Tokenize(";");
	    if(a){
		for(Int_t i = 0 ; i < a->GetEntries(); i++){
		    TObjString* s = (TObjString*) a->At(i);
		    TString t = s->GetString();
		    t.ReplaceAll(";","");

		    if( t.Contains("type" )       !=0){
			t.ReplaceAll("type=","");
			d->type = t;
		    } else if (t.Contains("channel"   )!=0){
			t.ReplaceAll("channel=","");
			d->channel = t;
		    } else if (t.Contains("T1"   )!=0 && t.Contains("fracT1"   )==0){
			t.ReplaceAll("T1=","");
			d->T1 = t.Atof();
		    } else if (t.Contains("T2"   )!=0){
			t.ReplaceAll("T2=","");
			d->T2 = t.Atof();
		    } else if (t.Contains("fracT1"   )!=0){
			t.ReplaceAll("fracT1=","");
			d->fracT1 = t.Atof();
		    } else if (t.Contains("blast")!=0){
			t.ReplaceAll("blast=","");
			d->blast = t.Atof();
		    } else if (t.Contains("A2"   )!=0){
			t.ReplaceAll("A2=","");
			d->A2 = t.Atof();
		    } else if (t.Contains("A4"   )!=0){
			t.ReplaceAll("A4=","");
			d->A4 = t.Atof();
		    } else if (t.Contains("V1"   )!=0){
			t.ReplaceAll("V1=","");
			d->V1 = t.Atof();
		    } else if (t.Contains("V2"   )!=0){
			t.ReplaceAll("V2=","");
			d->V2 = t.Atof();
		    } else if (t.Contains("eta"  )!=0){
			t.ReplaceAll("eta=","");
			d->eta = t.Atof();
		    }
		}
		delete a;

		//--------------------------------------------------------
		// add decay to the list odf decays for this energy
		it = datBase.mEtoDecay.find(energy);
		if(it == datBase.mEtoDecay.end()){ // create a new Set
		    vector<thermaldecay*> vd;
		    vd.push_back(d);
		    datBase.mEtoDecay [energy] = vd;
		} else {
		    it->second.push_back(d);
		}
		//--------------------------------------------------------
	    }
	}

	return kTRUE;
    }
};
// END DATABASE HELPERS
//----------------------------------------------------------------




Int_t pluto_embedded(TString outdir = "./", TString outfile = "test",TString type = "pi0",
		     Int_t nEvents = 1000, Int_t source = 1,Float_t Energy = 1.58, TString datbasefile = "database_mar19.dat",
		     TString vertexntuple = "")
{

    //-------------------------------------------------
    // RANDOM SEED
    // method;        0 (default) : /dev/random,
    // 1 : TRandom3 (local) with systime ,
    // 2 : TRandom3 (local) with systime + processid + iplast,
    // 3 : like 1 but gRandom,
    // 4 : like 2 but gRandom,
    // 5 : fixed (needs seed to be set)
    Int_t method = 0;
    HSeed hseed(method);
    Int_t seed = hseed.getSeed();

    PUtils::SetSeed(seed);
    UInt_t seed2 = (Int_t)PUtils::sampleFlat()*kMaxUInt;
    gRandom->SetSeed(seed2); // workarround for pluto BUG (new in 5.43 fixed in 5.45)
    //-------------------------------------------------

    //---------------CONIGURATION----------------------
    //

    Bool_t useVertex = kFALSE;
    if(vertexntuple.CompareTo("no") != 0 && vertexntuple.CompareTo("") != 0) useVertex = kTRUE;
    Int_t   sourceType     = source; // 0 : white single tracks , 1 : thermal source
    Int_t   asciiOut       = 0;    // write pluto ascci output for HGeant (==0 if we use HGeantOutput)
    Int_t   rootOut        = 0;    // write pluto root file
    Int_t   calcVertex     = useVertex ? 1 : 0;           // write the vertex to the ascii output for HGeant
    Bool_t  writeSeqNumber = useVertex ? kTRUE : kFALSE;  // write eventSeqNumber in addition
    Bool_t  writeIndex     = useVertex ? kTRUE : kFALSE;  // write parentIndex in addition

    //-------------------------------------------------
    // OUTPUT FILE
    outfile = gSystem->BaseName(outfile.Data());
    TString filename = Form("%s/%s",outdir.Data(),outfile.Data());
    filename.ReplaceAll("//","/");
    filename.ReplaceAll(".root","");
    filename.ReplaceAll(".evt","");
    //filename = Form("%s_pluto_%s",filename.Data(),type.Data());

    cout<<"############################################################"<<endl;
    cout<<"CONFIG OF PROGRAM:"<<endl;
    cout<<"output file   : "<<filename<<endl;
    cout<<"source switch : "<<sourceType<<endl;
    cout<<"calcvertex    : "<<calcVertex<<endl;
    cout<<"writeSeqNum   : "<<(Int_t)writeSeqNumber<<endl;
    cout<<"writeIndex    : "<<(Int_t)writeIndex<<endl;
    cout<<"############################################################"<<endl;
   //-------------------------------------------------

    PReaction* r = 0; // will be set later

    if(sourceType == 0)
    {   // white single tracks embedding

	cout<<"############################################################"<<endl;
	cout<<"CONFIG WHITE PARTICLE : "<<type<<endl;
	Float_t pmin               =   0.0;   // minimum mom  [GeV/c]
	Float_t pmax               = 2.000;   // maximum mom  [GeV/c]
	Float_t thetamin           =    10;   // minimum theta angle  [deg]
	Float_t thetamax           =    90;   // maximum theta angle  [deg]
	Float_t phimin             =     0;   // minimum phi angle    [deg]
	Float_t phimax             =    60;   // maximum phi angle    [deg]
	Float_t phistart           =     0;   // start phi angle (+ sector * 60. deg will be rotated later)
	Int_t numParticlePerSector =     1;   // number of Particles embedded per sector
        cout<<"pmin     = "<<setw(10)<<pmin    <<" pmax     = "<<setw(10)<<pmax<<endl;
	cout<<"thetamin = "<<setw(10)<<thetamin<<" thetamax = "<<setw(10)<<thetamax<<endl;
	cout<<"phimin   = "<<setw(10)<<phimin  <<" phimax   = "<<setw(10)<<phimax<<endl;
	cout<<"phistart = "<<setw(10)<<phistart<<" num Particle/sec = "<<numParticlePerSector<<endl;
	cout<<"############################################################"<<endl;

	PChannel** dummychannel  = 0;    // no channel needed here
	r = new PReaction(dummychannel,(Char_t*)filename.Data(),0,rootOut,0,calcVertex,asciiOut);
	
	PEmbeddedParticles* embedded = new PEmbeddedParticles();  //Construct the embedded container:

	PParticle* particle = new PParticle((Char_t*)type.Data(),1.,2.,3.); //Add an particle which we emit at a single point:
	embedded->AddParticle(particle); //Just add the particle to the container:

	if (1) { // nparticle / sector
	    embedded->SetSamplingSector(pmin, pmax, thetamin, thetamax, phimin, phimax, phistart, numParticlePerSector );
	} else { // nparticle / full acceptance
	    
	    embedded->SetSamplingSector(pmin, pmax, thetamin, thetamax,  -180.,   180., phistart, numParticlePerSector, 360., 1);
	}

	r->AddBulk(embedded);  //Add our container to the reaction:

    } else {  // thermal sources

	database datBase;
	if(!datBase.readDataBase(datBase,datbasefile)) return 1;
	//datBase.print();
        Bool_t ret = kFALSE;
	thermaldecay* d = datBase.getDecayParams(ret,Energy,type);
	if(d){
	    cout<<"############################################################"<<endl;
	    cout<<"SELECTED CHANNEL:"<<endl;
	    d->print();
	    cout<<"############################################################"<<endl;
	} else {
	    datBase.print();
	    exit(1);
	}
	if(d->type == "rho0"){ // code provided by S. Harabasz 28.11.2024
	    cout<<"############################################################"<<endl;
            cout<<"2 pion mass cut-off disabled for rho0"<<endl;
	    cout<<"############################################################"<<endl;
	    PFixedDecay *pmodel = new PFixedDecay("rho0_novdm_e-_e+", "Rho decay without VDM",-1);
	    pmodel->Add("rho0, parent");
	    pmodel->Add("e+, daughter");
	    pmodel->Add("e-, daughter");
	    makeDistributionManager()->Add(pmodel); //This enables the model
	    makeDistributionManager()->Enable("rho0_ee_e-_e+"); //Enable this for removing the pi-cutoff:
	    makeDistributionManager()->Print("decay_models");   //Finally check what we have done:
	}
        //makeDistributionManager()->Disable("eta_physics"); // Iza

	datBase.parseChannel(Energy,d);
	r = new PReaction(d->cfinal,(Char_t*)filename.Data(),d->vch.size(),rootOut,0,calcVertex,asciiOut);
    }

    r->SetWriteIndex(writeIndex);
    r->trackedParticles();
    r->Print();

    cout<<"############################################################"<<endl;
    cout<<"AFTER SETUP PREACTION"<<endl;
    cout<<"############################################################"<<endl;
    if(useVertex){
	cout<<"############################################################"<<endl;
	cout<<"USING VERTEX FILE : "<<vertexntuple<<endl;
	cout<<"############################################################"<<endl;
	PVertexFile *vertex = new PVertexFile(); //Construct the vertex container:
	vertex->OpenFile(vertexntuple);
	r->AddPrologueBulk(vertex); //add to prologue action
    } else {
	cout<<"############################################################"<<endl;
	cout<<"NO VERTEX SET "<<endl;
	cout<<"############################################################"<<endl;
    }

    PHGeantOutput* output = new PHGeantOutput();
    output->SetWriteSeqNumber(writeSeqNumber);
    output->OpenFile(Form("%s.evt",filename.Data()));
    r->AddFileOutput(output);

    r->setHGeant(0);

    if(useVertex){
	TFile *f = new TFile(vertexntuple.Data());
	if(f==NULL)
	{
	    cout<<"############################################################"<<endl;
	    cout << "file "<<vertexntuple.Data()<<" not found" << endl;
	    cout<<"############################################################"<<endl;
	    exit(1);

	}
	TNtuple *vertexnt = (TNtuple*)(f->Get("vertex"));
	if(vertexnt==NULL)
	{
	    cout<<"############################################################"<<endl;
	    cout << "NULL pointer to ntuple" << endl;
	    cout<<"############################################################"<<endl;
	    exit(1);
	}
	
	r->Loop(vertexnt->GetEntries()); // number of events
    } else {
	r->Loop(nEvents);
    }

    output->CloseFile();

    return 0;
}