#ifndef  ___PLUTO_DECAY__
#define  ___PLUTO_DECAY__

#include "PUtils.h"
#include "PParticle.h"
#include "PChannel.h"
#include "PReaction.h"


//----------------------------------------------------------------
// DATABASE HELPERS
typedef struct {
  Int_t   active;                 // if reaction is active
  Int_t   motherID;               // ID of the particle which should be decayed
  Double_t weight;                // weight for this decay (calculated by calcWeight())
  Double_t branching;             // branching ratio
  Int_t   sourceID;               // source id to uniquely identify the decay in the output
  TString type;                   // key word for this decay
  TString channel;                // channel description of this decay
  vector<PChannel*>  vch;         // list of channels
  vector<PParticle*> vm;          // mothers of decays
  vector<vector<PParticle*> > vd; // daughters of decays
  vector<PParticle*> vTracked;    // particles appearing in output
  PChannel** cfinal;              // the final PChannel of the PReaction
  PReaction* r; ;                 // PReaction pointer

  void clear(){
      active   = 0;
      motherID = -1;
      weight   = 1.;
      branching= 1.;
      sourceID = 0;
      type     = "";
      channel  = "";
      cfinal   = 0;
      vch.clear();
      vm.clear();
      for(UInt_t i = 0; i < vd.size(); i++) { vd[i].clear();}
      r = 0;
  }

  void print(Int_t i=-1, Int_t stat=-1){

      if(i < 0) {
	  cout<<setw(30)<<type
	      <<" chan "  <<setw(40)<<channel
	      <<" active : "<<active
	      <<" sourceID : "<<setw(12)<<sourceID
	      <<" weight : "<<setw(12)<<weight
	      <<" branch : "<<setw(12)<<branching<<flush;
              if(stat !=-1)cout<<",  number of decays = "<<stat<<flush;
	      cout<<endl;
      } else {
	  cout<<setw(3)<<i<<" : "<<setw(30)<<type
	      <<" chan "  <<setw(40)<<channel
	      <<" active : "<<active
	      <<" sourceID : "<<setw(12)<<sourceID
	      <<" weight : "<<setw(12)<<weight
	      <<" branch : "<<setw(12)<<branching<<flush;
              if(stat !=-1)cout<<",  number of decays = "<<stat<<flush;
	      cout<<endl;
      }
  }

} decay;

struct database {
    map<TString,decay*>         mDecay;          // map    decay keyword -> decay object
    map<Int_t, vector<decay*> > mIDDecayList;    // map    mother ID -> list of decays
    map<Int_t, vector<Int_t> > mIDDecayListStat; // map    mother ID -> list of stat of decays
    Int_t doRecoRho;                             // same as in params : setup for PLUTO part
    Int_t doRecoOmega;                           // doRecoXXX  re: reconstruct  mesons from smash to decay them to defined PLUTO channels
    Int_t doRecoPhi;
    Int_t doSecondDecay;
    Int_t doPluto;
    Int_t doTree;

    void print(){
        // print the content of the data base

	cout<<"#######################################################################################"<<endl;

	Int_t ct=1;
	cout<<"doRecoRho     : "<<doRecoRho<<endl;
	cout<<"doRecoOmega   : "<<doRecoOmega<<endl;
	cout<<"doRecoPhi     : "<<doRecoPhi<<endl;
	cout<<"doSecondDecay : "<<doSecondDecay<<endl;
	cout<<"doPluto       : "<<doPluto<<endl;
	cout<<"doTree        : "<<doTree<<endl;
	map<TString,decay*>::iterator it;
	for(it = mDecay.begin(); it != mDecay.end(); ++it ){
	    cout<<"---------------------------------------------------------------------------------------"<<endl;
	    it->second->print(ct);
	    ct++;
	}
	cout<<"---------------------------------------------------------------------------------------"<<endl;
	cout<<"#######################################################################################"<<endl;
    }

    void printStat(){
        // print the stat for the active channels

	cout<<"#######################################################################################"<<endl;
	cout<<"STAT OF PLUTO DECAYS:"<<endl;
	map<Int_t,vector<decay*> >::iterator it;
	map<Int_t,vector<Int_t> >::iterator it2;
	for(it = mIDDecayList.begin(); it != mIDDecayList.end(); ++it ){
	    cout<<"---------------------------------------------------------------------------------------"<<endl;
            it2 = mIDDecayListStat.find(it->first);
	    vector<decay*>& v = it->second;
	    vector<Int_t>& vstat = it2->second;
            for(UInt_t i = 0 ; i < v.size(); i++) {
		v[i]->print(i,vstat[i]);
	    }

	}
	cout<<"---------------------------------------------------------------------------------------"<<endl;
	cout<<"#######################################################################################"<<endl;
    }

    decay* getDecayParams(Bool_t& ret, TString type)
    {
	// return the decay object for a given type
	// from the database. Should be called after
        // readDataBase() has been called
	map<TString,decay*>::iterator it;

	it = mDecay.find(type);
	if(it != mDecay.end()){
	    ret = kTRUE;
	    return it->second;
	} else {
	    cout<<"ERROR : type = "<<type.Data()<<" not found in data base"<<endl;
	    ret = kFALSE;
	    return 0;
	}
    }

    Int_t getParticle(TString list,vector<PParticle*>& vp,Bool_t print=kFALSE)
    {
	// decomposes the comma separated list of PLUTO names
        // and creates PParticles which are returned in the vector vp
        // if print=kTRUE the names will be printed
	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(decay* d)
    {
	// takes the channel description from decay->channel
	// and decomposes it to PChannel setup and PParticles
	// for the given PReaction. All created obejcts are stored
        // in decay. PReactions will be created later.

	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
	d->motherID = d->vm [0]->ID();

	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-1;  // 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
	}
	//----------------------------------------------

	//----------------------------------------------
	// setting up the first level of decay chain connected
	// to the mother particle
	UInt_t nd = getParticle(l[0],vl[0],kFALSE);    // comma separated particle string -> PParticles
	d->vd.push_back(vl[0]);
	PParticle *moth = d->vm[0];                    // the mother of the decay is the main particle
	PParticle** smoth  = new PParticle* [nd+1];    // connect mother to Fireball
	smoth[0] = moth;
	for(UInt_t i = 0 ; i < nd; i++){               // connect the 1. level daughters
	    smoth[i+1] = vl[0][i];
	}
	PChannel* cmoth = new PChannel(smoth, nd, 1);  // create a channel for the first decay
	d->vch.push_back(cmoth);                       // add to the list of channels
	//----------------------------------------------


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

	    //----------------------------------------------
	    // analyze the level strings to build the decay channels
	    // for PReaction
	    for(Int_t i = 1 ; i < max+1; i++){
		getParticle(l[i],vl[i],kFALSE); // 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 = d->vm[i]->Name();
		    if( n != ml[i] ){
			cout<<"############################################################"<<endl;
			cout<<"ERROR : mismatch pid from index with stored name : pid from index = "<<d->vm[i]->Name()<<" 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
		//----------------------------------------------
	    }
	}

        //----------------------------------------------------
        // remember tracked particles
	for(UInt_t i = 0; i < d->vd.size(); i++){   // level
	    vector<PParticle*> vd = d->vd[i];
	    for(UInt_t j = 0 ; j < vd.size(); j++){ // loop daughter
		if(d->vm.size() >= i+1) {           // has at least one level more
                    Bool_t ismother = kFALSE;
		    for(UInt_t k = 0 ; k < d->vm.size(); k++){
			if( d->vm[k] == vd[j]) ismother = kTRUE;
		    }
                    if(!ismother) d->vTracked.push_back(vd[j]);

		} else { // final
                   d->vTracked.push_back(vd[j]);
		}
	    }
	}
        //----------------------------------------------------

        //----------------------------------------------------
	// creating the sourceID for this decay
	vector<Int_t> dID;
        for(UInt_t i = 0 ; i < d->vd[0].size() ; i++) dID.push_back(d->vd[0][i]->ID());
	d->sourceID = createSourceID(d->vm[0]->ID(),dID);
        //----------------------------------------------------

        //----------------------------------------------------
        // can this mother particle decay to different channels ?
	if(d->active){
	    map<Int_t, vector<decay*> >::iterator it;
	    it = mIDDecayList.find(d->motherID);
	    if(it == mIDDecayList.end()){  // new ID entry
		vector<decay*> v;
		vector<Int_t> vstat;
		v.push_back(d);
		vstat.push_back(0);
		mIDDecayList    [d->motherID] = v;
		mIDDecayListStat[d->motherID] = vstat;

	    } else { // add decay to list of decays for this ID
		vector<decay*>& v = mIDDecayList[d->motherID];
		v.push_back(d);
		vector<Int_t>& vstat = mIDDecayListStat[d->motherID];
		vstat.push_back(0);
	    }
	}
        //----------------------------------------------------


	//----------------------------------------------
	// build the final channel from collect channels above
	d->cfinal = new PChannel* [d->vch.size()];
	for(UInt_t i = 0 ; i < d->vch.size(); i++){
	    d->cfinal[i] = d->vch[i];
	}
	//----------------------------------------------

    }

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

		if( t.Contains("type" )       !=0){
		    t.ReplaceAll("type=","");
		    d->type = t;
                    ct++;
		} else if (t.Contains("channel"   )!=0){
		    t.ReplaceAll("channel=","");
		    d->channel = t;
                    ct++;
		} else if (t.Contains("active"    )!=0){
		    t.ReplaceAll("active=","");
		    d->active = t.Atoi();
                    ct++;
		} else if (t.Contains("branching" )!=0){
		    t.ReplaceAll("branching=","");
		    d->branching = t.Atof();
                    ct++;
		}
	    }
	    delete a;
	    if(ct >= 3) return kTRUE;
	}
	return kFALSE;
    }

    Bool_t readDataBase(database& datBase, parameters& params)
    {
	// read entries from predifined databasefile and
	// create a decay struct for each entry
	// format:
	// type=pi0dalitz ;channel=pi0:dilepton[e+,e-],g      ;active=1;
	// fileds are separated by ; and contain keyword=
	// channel description : use PLUTO names. , separate particles.
	// decays are described mother[daughter1,daughter2]. Decays
        // can be nested.
        // config keywords     keyword=value  [0 or 1]
	// defaults:
	//doRecoRho   = -1;
        //doRecoOmega = -1;
        //doRecoPhi   = -1;
        //doPluto     = -1;
        //doTree      = -1;
        // keywords in the databas file can be overwritten by command line parameters!
	//
	// empty lines, lines starting with '//' or '#' are skipped

	doRecoRho   = -1;
        doRecoOmega = -1;
        doRecoPhi   = -1;
        doSecondDecay = -1;
        doPluto     = -1;
        doTree      = -1;

	if(params.databasefile.EndsWith(".dat") == 1){

	    cout<<"############################################################"<<endl;
	    cout<<"READING DATABASE INPUT FILE : "<<params.databasefile<<endl;
	    cout<<"############################################################"<<endl;

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

	    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("doRecoRho")){
		    tmp.ReplaceAll("doRecoRho=","");
		    doRecoRho = tmp.Atoi();
		    continue;
		}

		if(tmp.Contains("doRecoOmega")){
		    tmp.ReplaceAll("doRecoOmega=","");
		    doRecoOmega = tmp.Atoi();
		    continue;
		}

		if(tmp.Contains("doRecoPhi")){
		    tmp.ReplaceAll("doRecoPhi=","");
		    doRecoPhi = tmp.Atoi();
		    continue;
		}

		if(tmp.Contains("doSecondDecay")){
		    tmp.ReplaceAll("doSecondDecay=","");
		    doSecondDecay = tmp.Atoi();
		    continue;
		}

		if(tmp.Contains("doPluto")){
		    tmp.ReplaceAll("doPluto=","");
		    doPluto = tmp.Atoi();
		    continue;
		}

		if(tmp.Contains("doTree")){
		    tmp.ReplaceAll("doTree=","");
		    doTree = tmp.Atoi();
		    continue;
		}

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

		if(decodeDecay(tmp,d)){
		    //--------------------------------------------------------
		    // add decay to the list odf decays for this energy
		    datBase.mDecay[d->type] = d;
		    parseChannel(d);
		    //--------------------------------------------------------
		} else {
		    cout<<"ERROR: readDatabase() Decay descriptor = "<<tmp<<" not a valid string!"<<endl;
		    exit(1);
		}
	    }

	} else{
	    cout<<"############################################################"<<endl;
	    cout<<"readDataBase() : data base input string "<<params.databasefile<<" not found!"<<endl;
	    cout<<"############################################################"<<endl;
	    TString  tmpall = params.databasefile;
	    TObjArray* a = tmpall.Tokenize("#");
	    if(a){
		for(Int_t i = 0 ; i < a->GetEntries(); i++){
		    TObjString* s = (TObjString*) a->At(i);
		    TString tmp = s->GetString();
                    tmp.ReplaceAll(" ","");

		    if(tmp.Contains("doRecoRho")){
			tmp.ReplaceAll("doRecoRho=","");
			doRecoRho = tmp.Atoi();
			continue;
		    }

		    if(tmp.Contains("doRecoOmega")){
			tmp.ReplaceAll("doRecoOmega=","");
			doRecoOmega = tmp.Atoi();
			continue;
		    }

		    if(tmp.Contains("doRecoPhi")){
			tmp.ReplaceAll("doRecoPhi=","");
			doRecoPhi = tmp.Atoi();
			continue;
		    }

		    if(tmp.Contains("doPluto")){
			tmp.ReplaceAll("doPluto=","");
			doPluto = tmp.Atoi();
			continue;
		    }

		    if(tmp.Contains("doSecondDecay")){
			tmp.ReplaceAll("doSecondDecay=","");
			doSecondDecay = tmp.Atoi();
			continue;
		    }

		    if(tmp.Contains("doTree")){
			tmp.ReplaceAll("doTree=","");
			doTree = tmp.Atoi();
			continue;
		    }

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

		    if(decodeDecay(tmp,d)){
			//--------------------------------------------------------
			// add decay to the list odf decays for this energy
			datBase.mDecay[d->type] = d;
			parseChannel(d);
			//--------------------------------------------------------
		    } else {
			cout<<"ERROR: readDatabase() Decay descriptor = "<<tmp<<" not a valid string!"<<endl;
			exit(1);
		    }

		}
	    }
	}

	if(params.doRecoRho == 0 && doRecoRho != -1) params.doRecoRho = doRecoRho;
	else if (params.doRecoRho == 1 &&  doRecoRho == 0 ) {
	    doRecoRho = params.doRecoRho;
	    cout<<"INFO: readDataBase() : option doRecoRho     has been overwritten by command line option!"<<endl;
	}

	if(params.doRecoOmega == 0 && doRecoOmega != -1) params.doRecoOmega = doRecoOmega;
	else if (params.doRecoOmega == 1 &&  doRecoOmega == 0 ) {
	    doRecoOmega = params.doRecoOmega;
	    cout<<"INFO: readDataBase() : option doRecoOmega   has been overwritten by command line option!"<<endl;
	}

	if(params.doRecoPhi == 0 && doRecoPhi != -1) params.doRecoPhi = doRecoPhi;
	else if (params.doRecoPhi == 1 &&  doRecoPhi == 0 ) {
	    doRecoPhi = params.doRecoPhi;
	    cout<<"INFO: readDataBase() : option doRecoPhi     has been overwritten by command line option!"<<endl;
	}

	if(params.doSecondDecay == 0 && doSecondDecay != -1) params.doSecondDecay = doSecondDecay;
	else if (params.doSecondDecay == 1 &&  doSecondDecay == 0 ) {
	    doSecondDecay = params.doSecondDecay;
	    cout<<"INFO: readDataBase() : option doSecondDecay has been overwritten by command line option!"<<endl;
	}

	if(params.doPluto == 0 && doPluto != -1) params.doPluto = doPluto;
	else if (params.doPluto == 1 &&  doPluto == 0 ) {
	    doPluto = params.doPluto;
	    cout<<"INFO: readDataBase() : option doPluto       has been overwritten by command line option!"<<endl;
	}

	if(params.doTree == 0 && doTree != -1) params.doTree = doTree;
	else if (params.doTree == 1 &&  doTree == 0 ) {
	    doTree = params.doTree;
	    cout<<"INFO: readDataBase() : option doTree        has been overwritten by command line option!"<<endl;
	}
	return kTRUE;
    }



    void calcWeights(){
	// If n decays for the same mother ID
	// exists each channel will be decayed with 1/n
        // the weight factor = n * branching will e stored in output
	map<Int_t, vector<decay*> >::iterator it;
	for(it = mIDDecayList.begin() ; it != mIDDecayList.end(); ++it){
	    vector<decay*>& v = it->second;
	    Double_t weight = v.size();
	    for(UInt_t i = 0 ; i < v.size(); i++){
		v[i]->weight = weight * v[i]->branching;
	    }
	}
    }


};
// END DATABASE HELPERS
//----------------------------------------------------------------


//----------------------------------------------------------------
//  Particle helpers
void printPParticle(PParticle* p)
{
    // print PParticle values for debug info
    cout<<"ID "<<setw(4)<<p->ID()<<" px "<<setw(12)<<p->Px()<<" py "<<setw(12)<<p->Py()<<" pz "<<setw(12)<<p->Pz()
	<<" M "<<setw(12)<<p->M()<<" E "<<setw(12)<<p->E()
	<<" genInfo "<<setw(12)<<p->GetSourceId()<<" genInfo1 "<<setw(12)<<p->GetParentId()<<" genInfo2 "<<setw(12)<<p->GetParentIndex()  <<endl;

}

Int_t smashToPParticle(smash_p& p,PParticle* pp)
{
    // fill PParticle from smash_p
    pp->SetPxPyPzE(p.px,p.py,p.pz,p.E);
    pp->SetVertex(p.x,p.y,p.z,p.form_time);
    pp->SetID(p.geantID);
    pp->SetSourceId(p.genInfo1);
    pp->SetParentId(p.genInfo2);
    pp->SetParentIndex(p.genInfo3);
    pp->SetW(p.weight);
    return 0;
}

Int_t PParticleToSmash(smash_p& p,PParticle* pp)
{
    // fill smash_p from PParticle
    p.reset();
    p.E         = pp->E();
    p.mass      = pp->M();
    p.p0        = pp->E();
    p.px        = pp->Px();
    p.py        = pp->Py();
    p.pz        = pp->Pz();
    p.pdg       = pp->ID();
    p.form_time = 100000;
    p.geantID   = pp->ID();
    p.charge    = pp->Charge();
    p.genInfo1  = pp->GetSourceId();
    p.genInfo2  = pp->GetParentId();
    p.genInfo3  = pp->GetParentIndex();
    p.weight    = pp->W();
    p.state     = 3;

    p.v.SetXYZM(p.px,p.py,p.pz,p.mass);

    mapGeant3ToPDG(obj.PDG,p);
    return 0;
}
//----------------------------------------------------------------


//----------------------------------------------------------------
// PLUTO Tree
void fillPutoTree(objects& ob, Double_t evtplane)
{
    // fill the event to a PLUTO tree (which has to be created before)
    obj.PLUTO_npart    = ob.evt.nParticle;
    obj.PLUTO_impact   = ob.evt.impact;
    obj.PLUTO_evtplane = evtplane;
    Int_t ct = 0;

    ob.PLUTO_evt->Clear();

    for(UInt_t i = 0 ; i < ob.particles.size(); i++)
    {
	smash_p& p = ob.particles[i];
	if(p.state <  0) continue;

	PParticle* pp = (PParticle*)ob.PLUTO_evt->ConstructedAt(ct);
	new(pp) PParticle();
        ct++;
	smashToPParticle(p,pp);
    }

    ob.PLUTO_tree->Fill();
}
//----------------------------------------------------------------

//----------------------------------------------------------------
//  decaying Particles with PLUTO
class InlineDecay {
public:
    UInt_t   seed;     // random seed
    database datBase;  // type -> decay

    InlineDecay(parameters& params)
    {
	PUtils::SetSeed(params.seed);

	makeDistributionManager()->Disable("eta_physics"); // according to Marten (Iza) , disto physics do messup

	makeDistributionManager()->Exec("vdm");
	makeDistributionManager()->Exec("elementary");
	makeDistributionManager()->Exec("brems: kaptari; weighting");
	makeDistributionManager()->Exec("brems: elastic");
	makeDistributionManager()->Exec("brems: fsi");
	makeDistributionManager()->Exec("dalitz_mod: krivoruchenko");
	//makeDistributionManager()->Exec("dalitz_mod: static_br_thresh=0.100 ; flat_generator"); // cause problem with w -> e+ e- pi0

	Float_t  hbar=6.582122e-25;// reduced Planck constant GeV*s
	Int_t pid_X17 = makeStaticData()->AddParticle(-1,"X17", 0.017);
	makeStaticData()->AddAlias("X17","X17");
	makeStaticData()->SetParticleTotalWidth("X17", hbar/4e-14);
	makeStaticData()->SetParticleParity("X17", 1);
	makeStaticData()->AddDecay("X17 --> e+ e-", "X17", "e+, e-", 1);
	makeStaticData()->AddDecay("eta --> pi+ pi- X17", "eta", "pi+, pi-, X17", 0.0001);


	makeStaticData()->AddDecay("eta --> e+ + e- + pi+ pi-", "eta", "e+,e-,pi+, pi-", 0.0001);

        datBase.readDataBase(datBase,params);

    }

    Bool_t isKnownDecay(Int_t id){
	// checks if there are setup decays for
        //the mother ID

       if (datBase.mIDDecayList.find(id) != datBase.mIDDecayList.end()) return kTRUE;
       else return kFALSE;
    }

    void initReactions()
    {
	// loops over data base and creates PReactions
        // for all active decays
	map<TString,decay*>::iterator it;
	for(it = datBase.mDecay.begin(); it != datBase.mDecay.end(); ++it){
	    decay* d = it->second;
            if(d->active != 1) continue;
 	    cout<<"############################################################"<<endl;
	    cout<<"------------------------------------------------------------"<<endl;
            cout<<"Building Reaction "<<d->channel<<endl;
 	    cout<<"------------------------------------------------------------"<<endl;
           TString dummy ="";
	    d->r = new PReaction(d->cfinal,(Char_t*)dummy.Data(),d->vch.size(),0/*rootOut*/,0, 0/*calcVertex*/,0/*asciiOut*/);
	    d->r->trackedParticles();
	    d->r->IsInline();
	    d->r->Print();
 	    cout<<"############################################################"<<endl;
	}
        datBase.calcWeights();
	datBase.print();

   };


   decay* decayPluto(objects& obj,Int_t nEvt, Int_t id, Double_t px,Double_t py, Double_t pz, Double_t E, vector<PParticle*>& daughters)
   {
       // calls PReaction for particle and returns a vector
       // of filled daughter PParticles
       daughters.clear();
       if(!isKnownDecay(id)) {
	   cout<<"ERROR in event "<<nEvt<<" no Decay found for id = "<<id<<endl;
	   return 0;
       }

       vector<decay*>& v = datBase.mIDDecayList[id];
       vector<Int_t>& vstat = datBase.mIDDecayListStat[id];
       UInt_t index = (UInt_t) (gRandom->Rndm() * v.size());
       if(index > (v.size() - 1)) index = v.size() - 1;
       decay* d = v[index];

       if(obj.debug > 0){
	   cout<<"------------------------------------------------------------"<<endl;
	   cout<<"decayPluto() ID "<<id<<" n decays "<<v.size()<<" index "<<index<<" decay "<<d->channel<<endl;
       }
       d->vm[0]->SetPxPyPzE(px,py,pz, E); // manipulate the source particle

       Int_t ret = d->r->Loop(1,0,0);     // decay particle with defined channel verbose off
       if(ret == 0) {
           d->vm[0]->SetPxPyPzE(px,py,pz, E);
	   cout<<"ERROR : decayPluto() : PLUTO decay failed! Skip this decay! ID "<<id<<" n decays "<<v.size()<<" index "<<index<<" decay "<<d->channel<<endl;

	   printPParticle(d->vm[0]);
	   return d;
       }
       vstat[index]++;

       for(UInt_t i = 0; i < d->vTracked.size(); i ++){ // copy daughters, but not source particle (vm[0])
           d->vTracked[i]->SetSourceId(d->sourceID);
           d->vTracked[i]->SetW(d->weight);
	   daughters.push_back(d->vTracked[i]);
       }

       return d;
   }



};


InlineDecay*  inlineDecay;
//----------------------------------------------------------------

#endif