#ifndef ___CONVERT_OSCAR13___
#define ___CONVERT_OSCAR13___


void printEvent(objects& obj,Int_t create,TString msg = " : ")
{
    vector<smash_p>& vpart = obj.particles;
    smash_p *p1;
    cout<<"-------------------------------------------------------------------------"<<endl;
    cout<<"Event "<<obj.evt.evtNr<<endl;
    cout<<msg<<" "<<obj.evt.nParticle<<" event full "<<vpart.size()<< " created "<< create<<endl;
    for (UInt_t i = 0; i < vpart.size() ; i++)
    {
	p1 = &vpart.at (i);
	cout<<setw(4)
	    <<i
	    <<" state "<<setw(2)    <<p1->state
	    <<" pdg "     <<setw(12)<<p1->pdg
	    <<" geantId " <<setw(6) <<p1->geantID
	    <<" charge "  <<setw(6) <<p1->charge
	    <<" ID "      <<setw(6) <<p1->ID
	    <<" time "    <<setw(12)<<p1->form_time
	    <<" time org "<<setw(12)<<p1->form_time_org
	    <<" moth1 "   <<setw(8) <<p1->pdg_mother1
	    <<" moth2 "   <<setw(8) <<p1->pdg_mother2
	    <<" info1 "   <<setw(12)<<p1->genInfo1
	    <<" info2 "   <<setw(12)<<p1->genInfo2
	    <<" info3 "   <<setw(12)<<p1->genInfo3
	    <<endl;

    }
}

Int_t createMother(objects& obj,Int_t srcID,smash_p& mother, vector<smash_p*> vd)
{

    mother.reset();

    for(UInt_t i = 0 ; i < vd.size(); i++ ) {
	if(vd[i]->state < 0) {
	    cout<<"ERROR: createMother() : daughter has negative state ! skip."<<endl;
            vd[i]->print(i);
	    mother.reset();
            printEvent(obj,0);
	    return -1;
	}
	mother.v += vd[i]->v;
    }
    mother.E      = mother.v.E();
    mother.mass   = mother.v.M();
    mother.p0     = mother.v.E();
    mother.px     = mother.v.Px();
    mother.py     = mother.v.Py();
    mother.pz     = mother.v.Pz();
    mother.pdg    = vd[0]->pdg_mother1;

    mother.geantID   = getPdgPIDtoGeant3(obj.PDG,mother.pdg,kFALSE);
    mother.charge    = (obj.PDG.GetParticle(mother.pdg)->Charge())/3;
    mother.form_time = vd[0]->form_time;
    mother.form_time_org = vd[0]->form_time_org;
    mother.weight    = 1.;
    mother.state     = 2;
    mother.genInfo1 = srcID;
    mother.v.SetXYZM(mother.px,mother.py,mother.pz,mother.mass);


    for(UInt_t i = 0 ; i < vd.size(); i++){
	vd[i]->state *= -1;
        vd[i]->form_time += 10000;
    }
    obj.evt.nParticle = obj.evt.nParticle - vd.size() + 1;

    return vd.size() - 1;
}

void  searchAndFlagDecaysReco(objects& obj, Bool_t print=kFALSE)
{
    // find decays and set geninfo respectivly
    // (only decays verified by mother ID + fomation time, not like omega from a before reconstructed rho)
    // Uses the decays defined in mDecay map

    

    vector<smash_p>& vpart = obj.particles;
    if(vpart.size() < 2) return;
    sort(vpart.begin(),vpart.end(),sortFormTime);
    vector<smash_p*> vd;
    vector<smash_p> vnew;
    vector<Int_t> vdaughters;
    vector<Int_t> vdaughtersID;
    UInt_t ich;
    Int_t srcID;

    Int_t oldN = vpart.size();

    smash_p *p1,*p2,*p3;
    Int_t i1,i2,i3;
    Int_t create = 0;
    Int_t nDau   = 0;
    Int_t nMoth  = 0;
    for (UInt_t i = 0; i < vpart.size()-1 ; i++)
    {
	p1 = &vpart.at (i);
	p2 = &vpart.at (i+1);
	i1 = 1;
        i2 = i + 1;
        i3 = -1;
	p3 = 0;
	if(i+2 < vpart.size()) p3 = &vpart.at (i+2);

	if(p1->pdg_mother1 != p2->pdg_mother1 || p1->form_time != p2->form_time ) continue;  // 2 daughter of decay
	vdaughters.clear();
	vd.clear();
	vdaughters.push_back(p1->pdg);
	vdaughters.push_back(p2->pdg);
        vd.push_back(p1);
        vd.push_back(p2);
	i++;

	if(p3){
	    if(p1->pdg_mother1 == p3->pdg_mother1 &&  p1->form_time == p3->form_time) {  // 3 daughters of decay
		vdaughters.push_back(p3->pdg);
		vd.push_back(p3);
                i3 = i2 + 1;
		i++;
	    }
	}

	if(obj.mDecay.findDecay(p1->pdg_mother1,ich,srcID,vdaughters))
	{
	    p1->genInfo1 = srcID;
	    p2->genInfo1 = srcID;
	    if(vdaughters.size() == 3) p3->genInfo1 = srcID;

	    if((params.doRecoRho   && isRho( p1->pdg_mother1)) ||
	       (params.doRecoOmega && p1->pdg_mother1 == 223)  ||
	       (params.doRecoPhi   && p1->pdg_mother1 == 333) )
	    {
		smash_p mother;
		createMother(obj,srcID,mother,vd);
		vnew.push_back(mother);
		create++;
                nDau += vdaughters.size();
                nMoth++;
	    }
	    //--------------------------------------------------

	    if(print)
	    {
		cout<<"-------------------------------------------------------------------------"<<endl;
		cout<<"Event "<<obj.evt.evtNr<<" mother ID "<< p1->pdg_mother1<<endl;
		obj.mDecay.printDecay(obj.PDG,p1->pdg_mother1,ich);

		cout<<setw(4)
		    <<i1
		    <<" d1 state "<<setw(2) <<p1->state
		    <<" pdg "     <<setw(12)<<p1->pdg
		    <<" geantId " <<setw(6) <<p1->geantID
		    <<" ID "      <<setw(6) <<p1->ID
		    <<" time "    <<setw(12)<<p1->form_time
		    <<" moth1 "   <<setw(8) <<p1->pdg_mother1
		    <<" moth2 "   <<setw(8) <<p1->pdg_mother2
		    <<" info1 "   <<setw(12)<<p1->genInfo1
		    <<" info2 "   <<setw(12)<<p1->genInfo2
		    <<" info3 "   <<setw(12)<<p1->genInfo3
		    <<endl;
		cout<<setw(4)
		    <<i2
		    <<" d2 state "<<setw(2) <<p2->state
		    <<" pdg "     <<setw(12)<<p2->pdg
		    <<" geantId " <<setw(6) <<p2->geantID
		    <<" ID "      <<setw(6) <<p2->ID
		    <<" time "    <<setw(12)<<p2->form_time
		    <<" moth1 "   <<setw(8) <<p2->pdg_mother1
		    <<" moth2 "   <<setw(8) <<p2->pdg_mother2
		    <<" info1 "   <<setw(12)<<p2->genInfo1
		    <<" info2 "   <<setw(12)<<p2->genInfo2
		    <<" info3 "   <<setw(12)<<p2->genInfo3
		    <<endl;
		if(vdaughters.size()==3 )
		    cout<<setw(4)
			<<i3
			<<" d3 state "<<setw(2) <<p3->state
			<<" pdg "     <<setw(12)<<p3->pdg
			<<" geantId " <<setw(6) <<p3->geantID
			<<" ID "      <<setw(6) <<p3->ID
			<<" time "    <<setw(12)<<p3->form_time
			<<" moth1 "   <<setw(8) <<p3->pdg_mother1
			<<" moth2 "   <<setw(8) <<p3->pdg_mother2
			<<" info1 "   <<setw(12)<<p3->genInfo1
			<<" info2 "   <<setw(12)<<p3->genInfo2
			<<" info3 "   <<setw(12)<<p3->genInfo3
			<<endl;
	    }
            //--------------------------------------------------
	} // decay found
    } // end loop

    for(UInt_t i = 0; i < vnew.size(); i++) vpart.push_back(vnew[i]);


    if(create > 0) { sort(vpart.begin(),vpart.end(),sortFormTime); }

    if(params.doRecoOmega)  // reco omega from  recontructed rho + pion
    {
	vnew.clear();
	create = 0;
        Int_t n = vpart.size();
	for (Int_t i = 0; i < n-1 ; i++)
	{
	    p1 = &vpart.at (i);
	    if(p1->state < 0) continue;
	    if(!(isRhoReco(p1->geantID) && p1->genInfo2 == 0)) continue; // reco rho 0,+,-

	    Int_t j = i;
	    if(j + 1 < n) j += 1;

	    for ( ; j >= 0; j--) { // starting from nearest neighbours, including adjacent 'later one' for safety.
		if ( j == i ) continue;
		if (j < 0) break;

		p2 = &vpart.at (j);

		if(p2->state < 0 || p1->state < 0) continue; // p1 might have changed state by createMother()

		if( p2->pdg_mother1 != 223 )      continue; // need daughter of omega
		if(!isPion (p2->pdg))             continue; // need pion
		if( p1->charge + p2->charge != 0) continue; // omega is neutral

		vector<smash_p*> vdtmp;
		vdtmp.push_back(p2);
		vdtmp.push_back(p1);

		vector<Int_t> vdID;
		vdID.push_back(p2->pdg);
		vdID.push_back(p1->pdg);

		smash_p omega;

		obj.mDecay.findDecay(p2->pdg_mother1,ich,srcID,vdID,kFALSE);
		createMother(obj,srcID,omega,vdtmp);

		create++;
                nDau += vdtmp.size();
                nMoth++;

		vnew.push_back(omega);
	    }
	}

	for(UInt_t i = 0; i < vnew.size(); i++) vpart.push_back(vnew[i]);
	if(create > 0) { sort(vpart.begin(),vpart.end(),sortFormTime); }
    } // doRecoOmega


    if(params.doRecoPhi)  // reco phi from  recontructed rho + pion
    {
	vnew.clear();
	create = 0;
        Int_t n = vpart.size();
	for (Int_t i = 0; i < n-1 ; i++)
	{
	    p1 = &vpart.at (i);
	    if(p1->state < 0) continue;
	    if(!(isRhoReco(p1->geantID) && p1->genInfo2 == 0)) continue;// reco rho 0,+,-

	    Int_t j = i;
	    if(j + 1 < n) j += 1;
	    for (; j >= 0; j--) { // starting from nearest neighbours, including adjacent 'later one' for safety.
		if ( j == i ) continue;
		if (j < 0) break;

		p2 = &vpart.at (j);

		if(p2->state < 0 || p1->state < 0) continue; // p1 might have changed state by createMother()

		if( p2->pdg_mother1 != 333 )      continue;
                if(!isPion(p2->pdg))              continue; // need pion
		if( p1->charge + p2->charge != 0) continue; // phi is neutral

		vector<smash_p*> vdtmp;
		vdtmp.push_back(p2);
		vdtmp.push_back(p1);

		vector<Int_t> vdID;
		vdID.push_back(p2->pdg);
		vdID.push_back(p1->pdg);

		smash_p phi;

		obj.mDecay.findDecay(p2->pdg_mother1,ich,srcID,vdID,kFALSE);
		createMother(obj,srcID,phi,vdtmp);

		create++;
                nDau += vdtmp.size();
                nMoth++;

		vnew.push_back(phi);
	    }
	}

	for(UInt_t i = 0; i < vnew.size(); i++) vpart.push_back(vnew[i]);
	if(create > 0) { sort(vpart.begin(),vpart.end(),sortFormTime); }
    } // doRecoPhi


    //---------------------------------------------------------
    // after reco rho than omega, some rho+- might be left usused
    // deactive the reco rho+- and put the original daughters of SMASH back
    if(params.doRecoRho)
    {
	Int_t reactivate = 0;
	vector<smash_p*> vrho;   // all rho+- state >0 in event
	vector<smash_p*> vrhoD;  // all daughter state <0 with rho+- mother
	for (UInt_t i = 0; i < vpart.size() ; i++)
	{
	    p1 = &vpart.at (i);
	    if(p1->state == 2 && isRhoReco(p1->geantID) && p1->genInfo2 == 0 && p1->charge != 0) { // reco rho+-, not used to reco omega
		vrho.push_back(p1);
	    }
	    if(p1->state < 0 && abs(p1->pdg_mother1) == 213) { // mother is 213 (rho+) or -213 (rho-)
		vrhoD.push_back(p1);
	    }
	}

	for (UInt_t i = 0; i < vrho.size() ; i++)
	{
	    p1 = vrho.at (i);

	    vector<smash_p*> vdfound;
	    for (UInt_t j = 0; j < vrhoD.size() ; j++) // check for each rho+- the corresponding daugthers
	    {
		smash_p* d = vrhoD[j];
		if(p1->form_time == d->form_time_org )  {
		    vdfound.push_back(d);
		}
	    }

	    if(vdfound.size() == 2) { // found both daughters of this rho+-

		for (UInt_t j = 0; j < vdfound.size() ; j++)
		{
		    vdfound[j]->state *= -1; // make them active again
		    vdfound[j]->form_time = vdfound[j]->form_time_org;
		}
		p1->state *=-1;      // deactivate rho+-
		obj.evt.nParticle ++;// correct nparticles
		reactivate ++;
	    }

	}

	if(reactivate > 0) { sort(vpart.begin(),vpart.end(),sortFormTime); }

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


    //---------------------------------------------------------
    // check constent number of active particles
    Int_t n = 0 ;
    for (UInt_t i = 0; i < vpart.size() ; i++)
    {
	p1 = &vpart.at (i);
	if(p1->state < 0) continue;
        n++;
    }

    if(obj.evt.nParticle != n) {
	cout<<"ERROR ::searchAndFlagDecaysReco() Number of particle inconsistent : evt.nParticle = "<<setw(4)<<obj.evt.nParticle
	    <<" actually : "<<setw(4)<<n
	    <<" old : "<<setw(4)<<oldN
	    <<" mother created : "<<setw(4)<<nMoth
	    <<" daughters removed : "<<setw(4)<<nDau
	    <<endl;
	TString msg = " wrong n particle : ";
	printEvent(obj,nMoth,msg);
    }
    //---------------------------------------------------------
}




Int_t convert_OSCAR13(objects& obj)
{
    // converts SMASH OSCAR13/OSCAR13 extended ascii format to HGeant ascii file.
    // Its assumed that the reaction is in lab frame already (see the smash settings)
    // K0_S and K0_L are equaly distributed 50:50
    // flag =  1;  // format Etot,px,py,pz,                ID,                              ,weight
    // flag = -2;  // format Etot,px,py,pz,                ID,sourceID,parentID,indexParent ,weight
    obj.evt.flag = -2; // output
    Int_t evtCt = 0;
    ifstream inFile(obj.input.Data());
    if(inFile.good())
    {
	smash_p   p;
        Char_t line[1000];


        //-----------------------------------------------------------
        // reading  smash header
	/*
	 #!OSCAR2013 particle_lists t x y z mass p0 px py pz pdg ID charge
	 # Units: fm fm fm fm GeV GeV GeV GeV GeV none none e
	 #
          or  extended format
         >= version 3.0
	 #!OSCAR2013Extended particle_lists t x y z mass p0 px py pz pdg ID charge ncoll form_time xsecfac proc_id_origin proc_type_origin time_last_coll pdg_mother1 pdg_mother2 baryon_number
	 # Units: fm fm fm fm GeV GeV GeV GeV GeV none none e none fm none none none fm none none none

	 < version 3.0
	 #!OSCAR2013Extended particle_lists t x y z mass p0 px py pz pdg ID charge ncoll form_time xsecfac proc_id_origin proc_type_origin time_last_coll pdg_mother1 pdg_mother2
         # Units: fm fm fm fm GeV GeV GeV GeV GeV none none e none fm none none none fm none none
	 #
	*/

	inFile.getline(line,1000);  // file format and vars decription
	TString l = line;
	if(!l.BeginsWith("#!OSCAR2013")) {
	    cerr<<"ERROR: input file is not OSCAR2013 format!"<<endl;
            return 1;
	}

	Bool_t extended = kFALSE;

	if(l.BeginsWith("#!OSCAR2013Extended")) {
	    extended = kTRUE;
	}

	inFile.getline(line,1000);  // units
	inFile.getline(line,1000);  // smash version
	TString l2 = line;
	l2.ReplaceAll("# SMASH-","");
        Float_t version = l2.Atof();
	cout<<"SMASH version : "<<version<<" !"<<endl;

        Int_t oldN = 0;
	//-----------------------------------------------------------

	while(!inFile.eof())
	{
	    //-----------------------------------------------------------
            // reading begin event
	    inFile.getline(line,1000);  // # event 1 out 435
	    if(!inFile.good()) break;
	    Int_t n = sscanf(line,"%*s %*s %d %*s %d",&obj.evt.evtNr,&obj.evt.nParticle);
	    //-----------------------------------------------------------
	    if(n==2){
		if(obj.debug >=3){
		    cout<<"-------------------------------------------------------------------------------"<<endl;
		    cout<<"Begin event :"<<setw(5)<<obj.evt.evtNr<<setw(4)<<" nParticle : "<<obj.evt.nParticle <<endl;
		}
                Bool_t ok=kTRUE;
		//-----------------------------------------------------------
		// particle eof event
		for(Int_t i = 1; i <= obj.evt.nParticle; i++ ){
                    p.reset();
		    if(!extended) {
			inFile.getline(line,1000);  //  t x y z mass p0 px py pz pdg ID charge
			n = sscanf(line,"%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d "
				   ,&p.t,&p.x,&p.y,&p.z,&p.mass,&p.p0,&p.px,&p.py,&p.pz,&p.pdg,&p.ID,&p.charge);

			if(n!=12){
			    cerr<<"ERROR: expected line is not a particle : "<<line<<endl;
			    ok = kFALSE;
			    break;
			}

		    } else {
			if(version >= 3){
			    inFile>> p.t >> p.x >> p.y >> p.z >> p.mass >> p.p0 >> p.px >> p.py >>p.pz >> p.pdg >> p.ID >> p.charge
				>> p.ncoll >> p.form_time >> p.xsecfac >> p.proc_id_origin >> p.proc_type_origin >> p.time_last_coll >> p.pdg_mother1 >> p.pdg_mother2 >> p.baryon_number;

			} else {
			    inFile>> p.t >> p.x >> p.y >> p.z >> p.mass >> p.p0 >> p.px >> p.py >>p.pz >> p.pdg >> p.ID >> p.charge
				>> p.ncoll >> p.form_time >> p.xsecfac >> p.proc_id_origin >> p.proc_type_origin >> p.time_last_coll >> p.pdg_mother1 >> p.pdg_mother2;
			}
			if(inFile.fail()){

			    cerr<<"ERROR: expected line is not a particle!"<<endl;
			    ok = kFALSE;
			    break;
			}
                        if(version >= 3 ) n = 21;
                        else              n = 20;
		    }

		    p.form_time_org = p.form_time;

		    mapPdgPIDtoGeant3(obj.PDG,p,kTRUE); // mass + charge + weight + geantID of p are set
                    setGenInfoSmash  (obj.PDG,p,kTRUE);
		    p.v.SetXYZM(p.px,p.py,p.pz,p.mass);

		    p.E = p.v.E();
		    p.state = 1;

		    //-----------------------------------------------------------
                    // do not store spectators
		    if(params.doIncludeSpectators == 0) {
			if( (extended && p.ncoll  == 0) || (!extended && p.t == 0.) ) {
			    p.state *= -1;
			} else obj.particles.push_back(p);
		    }
		    else obj.particles.push_back(p);
		    //-----------------------------------------------------------
		} // end particles
		//-----------------------------------------------------------

		if(!ok) {
		    break;
		}

		//-----------------------------------------------------------
		// reading end event
		inFile.getline(line,1000);  // example:  # event 1 end 0 impact   6.004
		TString t=line;
		if(t.Length()==0) inFile.getline(line,1000);  // example:  # event 1 end 0 impact   6.004
		Int_t ev;

		n = sscanf(line,"%*s %*s %d %*s %*d %*s %lf",&ev,&obj.evt.impact);
		//-----------------------------------------------------------

		if(n==2){
		    if(ev!=obj.evt.evtNr){
			cerr<<"ERROR: end event does not match begin event : begin = "<<obj.evt.evtNr<<" end = "<<ev<<endl;
                      //  break;
		    }
		    if(obj.debug >=3){
			cout<<"End  event :"<<setw(5)<<ev<<setw(4)<<" impact : "<<obj.evt.impact<<endl;
			cout<<"-------------------------------------------------------------------------------"<<endl;
		    }

		} else {
		    cerr<<"ERROR: expected line is not an end event : \""<<line<<"\""<<endl;
		    break;
		}

		obj.evt.nParticle = obj.particles.size(); // if spectators got removed

		//-----------------------------------------------------------
		// reconstruct Rho, Omega and Phi Mesons by their decay products
		// and replace decay products to later decay in PLUTO
		//
		// rho+-0-> 2 pions
		// omega -> pi0 gamma
		// omega -> rho pi
		// omega -> pi pi.
		// phi -> K+K-
		// phi -> K0 K0bar
		// phi -> eta gamma
		// phi -> pi0 gamma
		oldN = obj.particles.size();
                searchAndFlagDecaysReco(obj,kFALSE);
		//-----------------------------------------------------------

#ifdef WITHPLUTO
		//-----------------------------------------------------------
		// loop over event and decay mesons inline with PLUTO
                Bool_t create = 0;
                Int_t indexMoth = 0;
		if(params.doPluto == 1)
		{
		    vector<PParticle*>  daughters;
		    for(UInt_t i = 0 ; i < obj.particles.size(); i++)
		    {
			smash_p& p = obj.particles[i];
			Int_t id = p.geantID;
                        if(p.state <  0) continue;
			if(p.state == 3) continue;

			if(inlineDecay->isKnownDecay(id))
			{

			    if(obj.debug >= 2)
			    {
				cout<<"---------------------------------------------------------"<<endl;
				cout<<"PLUTO:  INPUT id "<<id<<" "<<p.p0<<" mass "<<p.mass<<endl;
				p.print(i);
			    }

			    decay* d = inlineDecay->decayPluto(obj,i, id , p.px, p.py, p.pz, p.p0,daughters);

			    if(d && daughters.size() > 0)
			    {
				if(obj.debug >= 2) cout<<"PLUTO: OUTPUT ndaughter "<<n<<" from "<<d->channel<<endl;
                                indexMoth++;
                                create++;
				p.state *= -1;
				UInt_t nd = daughters.size();
				for(UInt_t j = 0; j < nd; j ++)
				{
				    smash_p pd;
				    PParticle* dau = daughters[j];
				    dau->SetParentIndex(indexMoth);

				    PParticleToSmash(pd,dau);
				    obj.particles.push_back(pd);

				    obj.evt.nParticle ++ ;  // daughter add

				    if(obj.debug >= 2) pd.print(j);

				} // loop nd
				obj.evt.nParticle--; // mother remove
			    } // if d
			}
		    }

		    if(create > 0) sort(obj.particles.begin(),obj.particles.end(),sortFormTime);



		    //-----------------------------------------------------------
		    // looking for PLUTO decay output if a know decay channel
		    // needs to be decayed (example w -> e+ e- pi0    pi0 -> e+ e- g)
                    // the weight for this second generation decay has to calculated
		    if(params.doSecondDecay == 1)
		    {
			create = 0;
			daughters.clear();


			for(UInt_t i = 0 ; i < obj.particles.size(); i++)
			{
			    smash_p& p = obj.particles[i];
			    Int_t id = p.geantID;
			    if(p.state <  0) continue;
			    if(p.state != 3) continue;

			    if(inlineDecay->isKnownDecay(id))
			    {

				if(obj.debug >= 2)
				{
				    cout<<"---------------------------------------------------------"<<endl;
				    cout<<"SECOND DEACY PLUTO:  INPUT id "<<id<<" "<<p.p0<<" mass "<<p.mass<<endl;
				    p.print(i);
				}

				decay* d = inlineDecay->decayPluto(obj,i, id , p.px, p.py, p.pz, p.p0,daughters);

				if(d && daughters.size() > 0)
				{
				    if(obj.debug >= 2) cout<<"SECOND DECAY PLUTO: OUTPUT ndaughter "<<n<<" from "<<d->channel<<endl;
				    indexMoth++;
				    create++;
				    p.state *= -1;
				    UInt_t nd = daughters.size();
				    for(UInt_t j = 0; j < nd; j ++)
				    {
					smash_p pd;
					PParticle* dau = daughters[j];
					dau->SetParentIndex(indexMoth);

					PParticleToSmash(pd,dau);
					pd.weight *= d->weight;               // take into account weights
					pd.genInfo1 = p.genInfo1; // mark the daughters of second generation decay with source ID of the mother to allow reco later

					obj.particles.push_back(pd);

					obj.evt.nParticle ++ ;  // daughter add

					if(obj.debug >= 2) pd.print(j);

				    } // loop nd
				    obj.evt.nParticle--; // mother remove
				} // if d
			    }
			}

			if(create > 0) sort(obj.particles.begin(),obj.particles.end(),sortFormTime);
		    }
		    //-----------------------------------------------------------
		}
		//-----------------------------------------------------------

#endif
                if(obj.debug >=3) obj.printEvent("after event processing (no boost)");

		//-----------------------------------------------------------
		if(params.doBoost){
                   obj.boost();
		}
		//-----------------------------------------------------------

#ifdef WITHPLUTO
		//-----------------------------------------------------------
		// loop over event and fill PLUTO tree

		if(params.doTree == 1)
		{
                    fillPutoTree(obj,0.);
		}
		//-----------------------------------------------------------
#endif
		//-----------------------------------------------------------
                // check event size after reco + decay
		Int_t n = 0 ;
		for (UInt_t i = 0; i < obj.particles.size() ; i++)
		{
		    smash_p* p1 = &obj.particles.at (i);
		    if(p1->state < 0) continue;
		    n++;
		}

		if(obj.evt.nParticle != n) {
		    cout<<"ERROR ::convert_oscar13() Number of particle inconsistent : evt.nParticle = "<<setw(4)<<obj.evt.nParticle
			<<" actually : "<<setw(4)<<n
			<<" old : "<<setw(4)<<oldN
			<<endl;
		    TString msg = " wrong n particle : ";
		    printEvent(obj,-1,msg);
		}
		//-----------------------------------------------------------


		//-----------------------------------------------------------
                // stop if the requested number of events has been processed
		write_Header(obj);
		write_Event (obj);
		evtCt++;

		if(evtCt == obj.nEvts && obj.nEvts != -1 ) {
		    //goto STOP;
		    break;
		}
		//-----------------------------------------------------------

	    } else {
		cerr<<"ERROR: expected line is not an begin event : \""<<line<<"\""<<endl;
		break;
	    }
	}
    } else {
	cerr<<"ERROR: could not open : "<<obj.input<<"!"<<endl;
	return 1;
    }
//STOP:
#ifdef WITHPLUTO

    if(params.doRecoRho || params.doRecoOmega || params.doRecoPhi ) {
	inlineDecay->datBase.printStat();
    }

#endif
    inFile.close();
    obj.out.close();
    cout<<"Number of events read from file : "<<setw(6)<<evtCt<<" (requested "<<obj.nEvts<<")"<<endl;
    return 0;
}


#endif