#include "TString.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TRandom.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TTree.h"
#include "TVector3.h"
#include "TLorentzVector.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "TDecayChannel.h"


#include <fstream>
#include <iostream>
#include <iomanip>
#include <set>
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <vector>
#include <map>
#include <algorithm>


using namespace std;




#define WITHPLUTO        // needed for option --pluto (decay the reconstructed mesons inline using PLUTO)
#define VERSION 1.6
//     1.0 : initial commit to svn
//     1.1 : added K0L K0_SS handling  (50:50)
//     1.2 : added reconstruction of rho,omega and phi mesons from decays products for external decay (PLUTO)
//           added convert UrQmd 14 files
//           added  boost  CM -> LAB
//     1.3 : added decay of particle (PLUTO) , generic reaction build (leptons+hadrons)+ database input file
//           added geninfo, geninfo1, geninfo, weight to output (PLUTO to out put)
//     1.4 : in case recoRho etc, check event after omega reco for remaining rho+- and remove them(putting org smash particles in place)
//     1.5 : adding branching rations for weight calculations. added doSecondDecay option
//     1.6 : -incSpec option (by default spectators are not not propagated to the output)
//---------------------------------------------------------------------------------------------------------------
// helper structs + objects + conversion pid
#include "pid_converters.h"
#include "structs.h"
#include "pid_converters_struct.h"
//---------------------------------------------------------------------------------------------------------------


//---------------------------------------------------------------------------------------------------------------
// let pi0,eta,rho0,omega,phi inline decay with PLUTO
#ifdef WITHPLUTO
#include "pluto_decay.h"   // InlineDecay Class + pointer : requires PLUTO
#endif
//---------------------------------------------------------------------------------------------------------------


//---------------------------------------------------------------------------------------------------------------
// event output + converters for UrQMD (f14), smash (ROOT,OSCAR13)
#include "output.h"           // event output format
#include "convert_root.h"     // convert smash ROOT file input
#include "convert_urqmd.h"    // convert UrQMD f14 file input
#include "convert_oscar13.h"  // convert smash OSCAR13 input
//---------------------------------------------------------------------------------------------------------------



Int_t convert()
{

    if(gSystem->AccessPathName(obj.input.Data())){
	cerr<<"ERROR: input file : "<<obj.input<<" not found!"<<endl;
	return 1;
    }

    if(gSystem->AccessPathName(gSystem->DirName(obj.output.Data()))){
	cerr<<"ERROR: output dir : "<<gSystem->DirName(obj.output.Data())<<" not found!"<<endl;
	return 1;
    }

    obj.out.open(obj.output.Data());
    if(!obj.out.good()){
	cerr<<"ERROR: Could not open outfile : "<<obj.output<<"!"<<endl;
	return 1;
    }


    obj.out.precision(6);
    obj.out.setf( std::ios::fixed, std:: ios::floatfield );


    obj.PLUTO_out  = 0;
    obj.PLUTO_tree = 0;
    obj.PLUTO_evt  = 0;
    obj.PLUTO_npart= 0;
    obj.PLUTO_impact= 0;
    obj.PLUTO_evtplane= 0;

    if(params.doTree == 1)
    {

#ifdef WITHPLUTO

	TString name = obj.output;
        name.ReplaceAll(".evt","_pluto.root");


	obj.PLUTO_evt = new TClonesArray("PParticle", 1000);
	obj.PLUTO_out  = new TFile(name.Data(),"RECREATE");
	obj.PLUTO_tree = new TTree("data","data");
	obj.PLUTO_tree->Branch("Particles", &obj.PLUTO_evt     ,32000,99);
	obj.PLUTO_tree->Branch("Npart",     &obj.PLUTO_npart   , "Npart/I");
	obj.PLUTO_tree->Branch("Impact",    &obj.PLUTO_impact  , "Impact/D");
	obj.PLUTO_tree->Branch("Phi",       &obj.PLUTO_evtplane, "Phi/D");

#endif
    }




    if(obj.input.EndsWith(".root")){ // only availbale for smash
	cerr<<"input file is in smash root format!"<<endl;
	if(params.doRecoRho == 1) {
	    cerr<<"option doRecoRho for smash root file not suported! Abort."<<endl;
            return 1;
	}
	if(params.doRecoOmega == 1) {
	    cerr<<"option doRecoOmega for smash root file not suported! Abort."<<endl;
            return 1;
	}
	if(params.doRecoPhi == 1) {
	    cerr<<"option doRecoPhi for smash root file not suported! Abort."<<endl;
            return 1;
	}
	return convert_ROOT(obj);
    } else {
	ifstream inFile(obj.input.Data());
	if(inFile.good()){
	    Char_t line[1000];
	    inFile.getline(line,1000);  // file format and vars decription
	    TString l = line;
	    if(l.BeginsWith("#!OSCAR2013")) {
		cerr<<"input file is in OSCAR2013 format!"<<endl;

                Int_t ret = convert_OSCAR13(obj);
		if(params.doTree == 1)
		{
#ifdef WITHPLUTO
		    obj.PLUTO_tree ->Write();
                    obj.PLUTO_out -> Save();
		    obj.PLUTO_out -> Close();
#endif
		}

		obj.mDecay.printDecayMap(obj.PDG);

		return ret;
	    }
	    else if ( l.Contains("UQMD") && l.Contains("output_file  14"))
	    {
		cerr<<"input file is in UQMD f14 format!"<<endl;
		if(params.doRecoRho == 1) {
		    cerr<<"option doRecoRho for f14 file not suported! Abort."<<endl;
		    return 1;
		}
		if(params.doRecoOmega == 1) {
		    cerr<<"option doRecoOmega for f14 file not suported! Abort."<<endl;
		    return 1;
		}
		if(params.doRecoPhi == 1) {
		    cerr<<"option doRecoPhi for f14 file not suported! Abort."<<endl;
		    return 1;
		}
		return convert_UrQMD(obj);
	    } else {
                cerr<<"input file is in unknown format!"<<endl;
		return 1;
	    }
	}
    }

    return 0;
}

int main(int argc, char **argv)
{
    TROOT Analysis("Analysis","compiled analysis macro");

    //---------------------------------------------------------------------------------
    params.verbose_flag        = 0;
    params.nevents             = 100000000;
    params.help                = 0;
    params.doIncludeSpectators = 0;
    params.doRecoRho           = 0;
    params.doRecoOmega         = 0;
    params.doRecoPhi           = 0;
    params.doSecondDecay       = 0;
    params.doPluto             = 0;
    params.doTree              = 0;
    params.doBoost             = 0;
    params.seed                = 0;
    params.databasefile        ="";
    params.infile              ="";
    while (1)
    {
	static struct option long_options[] =
	{
	    // These options set a flag.
	    {"verbose"    ,  required_argument, 0, 'v'},
	    {"help"       ,  no_argument      , 0, 'h'},
	    {"input"      ,  required_argument, 0, 'i'},
	    {"output"     ,  required_argument, 0, 'o'},
	    {"events"     ,  required_argument, 0, 'e'},
	    {"incSpec"    ,  no_argument      , 0, 'a'},
	    {"recoRho"    ,  no_argument      , 0, 'r'},
	    {"recoOmega"  ,  no_argument      , 0, 'w'},
	    {"recoPhi"    ,  no_argument      , 0, 'p'},
	    {"ntuple"     ,  no_argument      , 0, 'n'},
	    {"secondDecay",  no_argument      , 0, 'c'},
	    {"pluto"      ,  no_argument      , 0, 'u'},
	    {"boost"      ,  required_argument, 0, 'b'},
	    {"seed"       ,  required_argument, 0, 's'},
	    {"databasefile", required_argument, 0, 'd'},
	    {0, 0, 0, 0}
	};
	// getopt_long stores the option index here.
	Int_t option_index = 0;

	Int_t c = getopt_long (argc, argv, "i:o:e:b:s:v:d:harwcpun", long_options, &option_index);

	// Detect the end of the options.
	if (c == -1)
	    break;

	switch (c)
	{
	case 0:
	    // If this option set a flag, do nothing else now.
	    if (long_options[option_index].flag != 0)
		break;
	    printf ("option %s", long_options[option_index].name);
	    if (optarg) printf (" with arg %s", optarg);
	    printf ("\n");
	    break;

	case 'v':
	    params.verbose_flag = atoi(optarg);
	    break;

	case 'a':
	    params.doIncludeSpectators = 1;
	    break;

	case 'r':
	    params.doRecoRho = 1;
	    break;

	case 'w':
	    params.doRecoOmega = 1;
	    break;

	case 'p':
	    params.doRecoPhi = 1;
	    break;

	case 'c':
	    params.doSecondDecay = 1;
	    break;

	case 'u':
	    params.doPluto = 1;
	    break;

	case 'n':
	    params.doTree  = 1;
	    break;

	case 'b':
	    params.doBoost = 1;
	    params.boostConfig = optarg ;
	    break;

	case 's':
	    params.seed = atoi(optarg);
	    break;

	case 'h':
	    params.help = 1;
	    break;

	case 'i':
	    params.infile = optarg ;
	    break;

	case 'o':
	    params.outfile = optarg ;
	    break;

	case 'e':
	    params.nevents = atoi(optarg) ;
	    break;

	case 'd':
	    params.databasefile = optarg ;
	    break;

	case '?':
	    // getopt_long already printed an error message.
	    break;

	default:
	    abort ();
	}
    }


    if(params.doPluto == 1 && params.databasefile == ""){
	cerr<<"ERROR: option --pluto requires  a database file"<<endl;
        return 1;
    }

    if(params.doPluto == 0 && params.doRecoRho == 1){
	cerr<<"ERROR: option --recoRho requires --pluto "<<endl;
        return 1;
    }

    if(params.doPluto == 0 && params.doRecoOmega == 1){
	cerr<<"ERROR: option --recoOmega requires --pluto "<<endl;
        return 1;
    }

    if(params.doPluto == 0 && params.doRecoPhi == 1){
	cerr<<"ERROR: option --recoPhi requires --pluto "<<endl;
        return 1;
    }

    if(params.doPluto == 0 && params.doSecondDecay == 1){
	cerr<<"ERROR: option --secondDecay requires --pluto "<<endl;
        return 1;
    }



    // Print any remaining command line arguments (not options).
    if (optind < argc)
    {
	cerr<<"ERROR: non-option ARGV-elements: ";
	while (optind < argc)
	    cerr <<" "<< argv[optind++];
	cerr<<endl;
	return 1;
    }
    //---------------------------------------------------------------------------------



    if(!params.help && params.infile != "" && params.outfile != ""){
	params.print();

	//--------------------------------------------------------
	// get real random seed (non blocking) or fixed seed from input option
	if(params.seed == 0){
	    UInt_t seed = 0;
	    Int_t fFileHandle = open("/dev/urandom", O_RDONLY);
	    Int_t ret = read(fFileHandle, &seed, sizeof(seed));
	    if(ret < 0){
		cerr<<"ERROR: Could not read from /dev/urandom!"<<endl;
		return 1;
	    }
	    else {
		cout<<"Random seed = "<<seed<<endl;
	    }
	    close(fFileHandle);
            params.seed = seed;
	    gRandom->SetSeed(seed);
	} else {
	    gRandom->SetSeed(params.seed);
	}
	//--------------------------------------------------------


#ifdef  WITHPLUTO
	    if(params.doPluto == 1 || params.databasefile !="") {
		inlineDecay = new InlineDecay(params);
		inlineDecay ->initReactions();
	    } else  inlineDecay = 0;
#endif

	obj.evt.eBeam = -1.0;   // for the moment
	obj.debug  = params.verbose_flag;
	obj.input  = params.infile;
	obj.output = params.outfile;
	obj.nEvts  = params.nevents;
	obj.init();

	//---------------------------------------------------------------------------------
	// analyze the boostOption string
	if(params.doBoost == 1)
	{
	    obj.evt.eBeam = params.boostConfig.Atof();
	    obj.evt.beta  = obj.evt.calcBeta(obj.evt.eBeam);
	    obj.evt.gamma = obj.evt.calcGamma(obj.evt.beta);

	    cout<<"----------------------------------------------"<<endl;
	    cout<<"Events will be boosted from CM -> LAB"<<endl;
	    cout<<"Ebeam = "<<obj.evt.eBeam<<"/u [GeV] "<<endl;
	    cout<<"beta  = "<<obj.evt.beta<<endl;
	    cout<<"gamma = "<<obj.evt.gamma<<endl;
	    cout<<"----------------------------------------------"<<endl;

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

	return convert();
    } else {
	cerr<<"version "<<VERSION<<endl;
	cerr<<"USAGE:  "<<argv[0]<<" -i infile -o outfile [-h,v,u,a,r,w,c,p,n] [-e nevents  num] [-b beam_energy ] [-s seed] [-d databasefile]"<<endl;
	cerr<<"-i (--input) infile  : smash OSCAR13/OSCAR13Extended ascii/Root file format or UrQmd f14 input (Center of mass system)."<<endl;
	cerr<<"-o (--output) outfile: will be in HGeant ascii evt format."<<endl;
	cerr<<"-v (--verbose)       : verbose output level (0 default, 1,2,3)"<<endl;
	cerr<<"-e (--events) num    : run num events (default is all events)"<<endl;
	cerr<<"-h (--help)          : print this help"<<endl;
	cerr<<"-a (--incSpec)       : include spectators (no default)"<<endl;
	cerr<<"-r (--recoRho)       : reconstruct rho mesons for external decay (no default, needs OSCAR13 format)"<<endl;
	cerr<<"-w (--recoOmega)     : reconstruct rho,omega mesons for external decay (no default, needs OSCAR13 format)"<<endl;
	cerr<<"-p (--recoPhi)       : reconstruct phi mesons for external decay (no default, needs OSCAR13 format)"<<endl;
	cerr<<"-c (--secondDecay)   : if PLUTO decays are used, known channels daughters will be be decayed too (no default, needs OSCAR13 format, compile #define WITHPLUTO )"<<endl;
	cerr<<"-u (--pluto)         : required my reco options (no default, needs OSCAR13 format, compile #define WITHPLUTO )"<<endl;
	cerr<<"-n (--ntuple)        : produce PLUTO Tree (no default)"<<endl;
	cerr<<"-b (--boost)         : boost CM to LAB frame (no default) format : \"EKin/u\" [GeV]"<<endl;
	cerr<<"-s (--seed)   num    : run with this random seed (no default)"<<endl;
	cerr<<"-d (--databasefile)  : decay data base file (required when -p is used)"<<endl;
        cerr<<endl;
        cerr<<endl;
	cerr<<"HINT: options recoRho,recoOmega,recoPhi,pluto,ntuple can be set in databasefile in the form keyword=value (0 or 1)"<<endl;
        cerr<<"      options set in the databasefile will be overwritten by command line options."<<endl;
	cerr<<endl;
	cerr<<endl;
	cerr<<"FORMAT database file example lepton decays using the reco options:"<<endl;
        cerr<<"Extension has to be .dat (active=0 decays will be ignored)"<<endl;
	cerr<<"#------------------------------------------------------------------"<<endl;
	cerr<<"doRecoRho=1"<<endl;
	cerr<<"doRecoOmega=1"<<endl;
	cerr<<"doRecoPhi=1"<<endl;
	cerr<<"doPluto=1"<<endl;
	cerr<<"doSecondDecay=1"<<endl;
	cerr<<"doTree=0"<<endl;
	cerr<<"dtype=pi0dalitz ;channel=pi0:dilepton[e+,e-],g;branching=1.198e-2  ;active=1;"<<endl;
	cerr<<"dtype=etadalitz ;channel=eta:dilepton[e+,e-],g;branching=0.6e-2    ;active=1;"<<endl;
	cerr<<"dtype=eta       ;channel=eta:e+,e-                                 ;active=0;"<<endl;
	cerr<<"dtype=rho0      ;channel=rho0:e+,e-           ;branching=0.00448e-2;active=1;"<<endl;
	cerr<<"dtype=w         ;channel=w:e+,e-              ;branching=0.00715e-2;active=1;"<<endl;
	cerr<<"dtype=wdalitz   ;channel=w:dilepton[e+,e-],pi0;branching=0.059e-2  ;active=1;"<<endl;
	cerr<<"dtype=phi       ;channel=phi:e+,e-            ;branching=0.0299e-2 ;active=1;"<<endl;
	cerr<<"#------------------------------------------------------------------"<<endl;
	cerr<<endl;
	cerr<<endl;
	cerr<<"FORMAT database as string example (line seperator \"#\"):"<<endl;
	cerr<<"option -d \"doRecoRho=1#doRecoOmega=1#doRecoPhi=1#doPluto=1#doTree=0#type=pi0dalitz;channel=pi0:dilepton[e+,e-],g;active=1;\""<<endl;
        cerr<<"The key word \"branching=\" for branching ratio is optional, by default branching is assumed to be 1."<<endl;
        cerr<<"The branching ratio is used when calculating the weights of a particle :  weight = active nchannels for this mother particle * branching."<<endl;
        cerr<<"For the above example w+wdalitz = 2 active channels. weight for daughters of w = 2*0.00715e-2."<<endl;
	cerr<<endl;
	cerr<<endl;
	cerr<<"HELP on geninfo, geninfo1, geninfo2 in output to identify sources:"<<endl;
	cerr<<"Standard smash output:"<<endl;
	cerr<<"   1. single particles keep"<<endl;
        cerr<<"      Int_t sourceID = pow(10,1*2) * parentID1"<<endl;
	cerr<<"      sourceID      += pow(10,0*2) * parentID2"<<endl;
	cerr<<"      genInfo1 = sourceID"<<endl;
	cerr<<"      genInfo2 = 1000000 + proc_type_origin"<<endl;
	cerr<<"      genInfo3 = 1000000 + proc_id_origin"<<endl;
	cerr<<"      example proton from elastic scattering :      1414 100001 100001"<<endl;
	cerr<<"   2. decay daughters from smash rho (0,+,-), w, phi (without reco options)"<<endl;
	cerr<<"      sourceID = motherID * pow(10,1*3) + daughter1ID * pow(10,1*2) + daughter2ID * pow(10,1*1) "<<endl;
 	cerr<<"      genInfo2 = 1000000 + proc_type_origin"<<endl;
	cerr<<"      genInfo3 = 1000000 + proc_id_origin"<<endl;
	cerr<<"      example :  rho+ -> pi+ p-    420807 100005 100003"<<endl;
	cerr<<"   3. decay daughters from PLUTO decays"<<endl;
	cerr<<"      sourceID = motherID * pow(10,1*3) + daughter1ID * pow(10,1*2) + daughter2ID * pow(10,1*1) "<<endl;
 	cerr<<"      genInfo2 = like PLUTO"<<endl;
	cerr<<"      genInfo3 = counter for each PLUTO decay in event "<<endl;
	cerr<<"      example :  pi0 -> dipleton[e+,e-],g    75101 7501(e+,e-) or 7(g) 1"<<endl;
	cerr<<endl;
	cerr<<endl;
	cerr<<"OUTPUT FORMAT from OSCAR13 extended input contains generator infos"<<endl;
	cerr<<"Event Header  : Evt Nr  NParticle Collision energy  Impact parameter HGeant format flag "<<endl;
	cerr<<"Particle line : E, px, py, pz, geantID, genInfo1, genInfo2, genInfo3, weight"<<endl;
	cerr<<"5	 6	 4.500000	 0.396000	 -2 "<<endl;
	cerr<<"1.311937     0.505859     0.246233     0.724438           13         1313      1000003      1000001     0.000000"<<endl;
	cerr<<"2.010489    -0.408474     0.391084     1.924376            8          808      1000005      1000002     0.000000"<<endl;
	cerr<<"2.473349     0.123280    -0.655407     2.189259           14         1414      1000005      1000003     0.000000"<<endl;
	cerr<<"0.132957     0.003891    -0.009491     0.132561            1        75101            7            1     0.011980"<<endl;
	cerr<<"0.027938    -0.013071     0.001796     0.024621            2        75101         7051            1     0.011980"<<endl;
	cerr<<"0.417129    -0.211485     0.025785     0.358616            3        75101         7051            1     0.011980"<<endl;
	cerr<<endl;
	cerr<<endl;
	cerr<<"OUTPUT FORMAT from UrQMD f14 input contains no generator infos"<<endl;
	cerr<<"Event Header  : Evt Nr  NParticle Collision energy  Impact parameter HGeant format flag "<<endl;
	cerr<<"Particle line : E, px, py, pz, geantID, weight"<<endl;
	cerr<<endl;
	cerr<<endl;

	params.print();

	return 1;
    }

}