#include "Generate.h"
#include "Setup.H"
#include "TDecayMode.H"
#include "HEPParticle.H"
#include "HEPEvent.H"
#include "HEPEVTEvent.H"
#include "PYJETSEvent.H"
#include "LUJETSEvent.H"
#include "HerwigEvent.H"
#include "GenerationDescription.H"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <unistd.h>
#include <assert.h>
#include "MCTesterEvent.H"
#include <THEPParticle.H>

#include "UserEventAnalysis.H"
#include "UserTreeAnalysis.H"
// ROOT-specific include files

#include "TFile.h"
#include "TInterpreter.h"
#include "TMethodCall.h"
#include "TSystem.h"

// global variables
TFile *file=0;

char decaymodes[2000][100];
int n_decaymodes=0;


GenerationDescription *setup_copy=0;


TMethodCall* userTreeAnalysis=0;


void MC_Initialize()
{
    
    

    if (!gInterpreter) {
	printf("FATAL ERROR: root not initialized\n");
	exit(-1);
    }

    if (Setup::stage==0) Setup::stage=1;

    // now try to open SETUP file.
    FILE *f=fopen("SETUP.C","r");
    if (f) {
	fclose(f);
	gInterpreter->ExecuteMacro("SETUP.C");
    } else {
	printf("\n\n");
	printf("	#############################\n");
	printf("	#                           #\n");
	printf("	# WARNING! NO SETUP.C file. #\n");
	printf("	#                           #\n");
	printf("	#############################\n");
	printf("\n");
    }

    printf("\n");
    printf(" *************************************\n");
    printf(" *              MC-TESTER            *\n");
    printf(" *-----------------------------------*\n");
    printf(" *    Testing decays of:%.8s         *\n",HEPParticle::GetParticleName(Setup::decay_particle));
    printf(" *                                   *\n");
    printf(" *   (c) Piotr    Golonka,   (1,2,3) *\n");
    printf(" *       Tomasz   Pierzchala,(2)     *\n");
    printf(" *       Zbigniew Was        (2,3)   *\n");
    printf(" *                                   *\n");
    printf(" *  1) FNPT,UMM, Krakow, Poland      *\n");
    printf(" *  2)      INP, Krakow, Poland      *\n");
    printf(" *  3)     CERN, Geneva, Switzerland *\n");
    printf(" *************************************\n\n");


    if (Setup::EVENT == 0) {
	Setup::EVENT=&HEPEVT;
	printf(" -> Using default event record.\n");
    }
    if (Setup::EVENT== &HEPEVT) 
	printf(" -> Event record format: HEPEVT\n");
    else if (Setup::EVENT== &LUJETS) 
	printf(" -> Event record format: LUJETS\n");
    else if (Setup::EVENT== &PYJETS) 
	printf(" -> Event record format: PYJETS\n");
    else if (Setup::EVENT== &HerwigEVT) 
	printf(" -> Event record format: HerwigEVT\n");
    else if (Setup::EVENT== &MCTEVT) 
	printf(" -> Event record format: MCTEVT\n");
    /** else if (Setup::EVENT== &HepMCEvt) 
	printf(" -> Event record format: HepMCEvt\n");**/
    else
	printf(" -> Event record format: unknown\n");
	
    if (Setup::user_event_analysis) {
	printf(" -> Using User Event Analysis Code from object:%s\n",
			Setup::user_event_analysis->GetName());
    }


    if (Setup::UserTreeAnalysis) {
	// setup user tree analysis...
	printf(" -> User Tree Analysis routine:%s\n",Setup::UserTreeAnalysis);
	userTreeAnalysis=new TMethodCall(Setup::UserTreeAnalysis,"Setup::UTA_particle,Setup::UTA_partlist,Setup::UTA_nparams, Setup::UTA_params");
    }



    char fname[2048];

    if (Setup::stage==1) {
	if (Setup::result1_path==0 || strlen(Setup::result1_path)==0)
	    Setup::result1_path="mc-tester.root";

	Setup::ResolvePath(Setup::result1_path,fname);
        printf(" -> results from stage1 goes to:\n     %s\n",fname);

    } else {
	if (Setup::result2_path==0 || strlen(Setup::result2_path)==0)
	    Setup::result2_path="mc-tester.root";

	Setup::ResolvePath(Setup::result2_path,fname);
        printf(" -> results from stage2 goes to:\n     %s\n",fname);	

    }
    
    printf("\n");
	setup_copy=new GenerationDescription(Setup::setup);
	getcwd(setup_copy->gen_path,128);


    if (Setup::stage==2) {
	sprintf(setup_copy->gen_desc_1,"%s",Setup::gen2_desc_1);
	sprintf(setup_copy->gen_desc_2,"%s",Setup::gen2_desc_2);
	sprintf(setup_copy->gen_desc_3,"%s",Setup::gen2_desc_3);
    } else {
	sprintf(setup_copy->gen_desc_1,"%s",Setup::gen1_desc_1);
	sprintf(setup_copy->gen_desc_2,"%s",Setup::gen1_desc_2);
	sprintf(setup_copy->gen_desc_3,"%s",Setup::gen1_desc_3);
    }

    file=TFile::Open(fname,"RECREATE");

    if (!file) {
	printf(" ! ERROR: cannot open output file: %s\n",fname);
	exit(-1);
    }
    
    file->cd();


}


void MC_Finalize()
{
    // append user-mode histograms if any...
    if (Setup::user_histograms->GetEntries()) {
	file->cd();
	TDirectory *dir=file->mkdir("USER_HISTOGRAMS");
	dir->cd();
	dir->Append(Setup::user_histograms);
    }
    file->cd();


    file->Append(setup_copy);
    file->Write();

    printf("-------------END OF MC-TESTER RUN-------------------\n");

    int nchannels=TDecayMode::DecayModes->GetEntriesFast();
    printf("Total: %i channels found, %li events analyzed\n",nchannels, Setup::events_cnt);
    
    TIter next(TDecayMode::DecayModes);
    
    while (TDecayMode *dm=(TDecayMode*)next()) {
    
	printf("  %s  (%li entries)\n",dm->GetName(),dm->GetNEntries());
    
    }
    printf("\n Total entries: %li\n\n",TDecayMode::NFills);
}

void MC_Analyze()
{
  MC_Analyze(Setup::decay_particle);
}

void MC_Analyze(HEPEvent * event, double weight)
{
  MC_Analyze(Setup::decay_particle, weight, event);
}

void MC_Analyze(int particle_PDG, double weight, HEPEvent * event)
{

  if(particle_PDG!=Setup::decay_particle){ //print error message is 
    std::cout << "MC-TESTER ERROR: Decay particle pdg id has change since initialization. Please check that" << std::endl;
    std::cout << "Setup::decay_particle is set in SETUP.C and is consistent with the pdg id given to MC_Analyze()" << std::endl;
    exit(-1);
  }

#define MCDEBUG(statement) if (Setup::debug_mode) {std::cout<<"DEBUG: Event "<<EVENT->GetEventNumber()<<" ["<<	Setup::events_cnt <<"]: "<<statement<<std::endl;}
#define MCDEBUGPART(statement,particle) if (Setup::debug_mode) {std::cout<<"DEBUG: Event "<<EVENT->GetEventNumber()<<" ["<<	Setup::events_cnt <<"] :"<<statement<<std::endl; std::cout<<"      ";particle->ls();}


    int ndecayproducts=0;
/*
    int myParticles[10000];
    int nparticles=0;
    int i=0;
    int children[10000];
    int nchildren=0;
    int decayproducts[1000];
    HEPParticle *p=0;
*/

    if(event)
      Setup::EVENT=event;

    HEPEvent *EVENT=Setup::EVENT;
    HEPEvent *EVENT_ORIG=Setup::EVENT;

    if (Setup::user_event_analysis) {
	UserEventAnalysis &a=*(Setup::user_event_analysis);
	a.SaveOriginalEvent(EVENT);
	EVENT=a.ModifyEvent(EVENT);
	Setup::EVENT=EVENT;
    }

    //    std::cout << "MC-TESTER sees the event:"<<std::endl;
    //    EVENT->ls();

    Setup::events_cnt++;

    // firstly: find all particles of that kind
    HEPParticleList plist;
    EVENT->FindParticle(particle_PDG,&plist);
    if (plist.empty()){
	
	MCDEBUG("no "<< HEPParticle::GetParticleName(particle_PDG)<<" particle found");
	return;
    } else {
	MCDEBUG(plist.size()<<" particles found");
    }
/*
    for (i=1; i<=EVENT->GetNumOfParticles();i++) {
	p=EVENT->GetParticle(i);
	// Check if it decays:
	if (! p->Decays() ) continue;
	if ( p->GetPDGId() == particle_PDG ) {
	    // protection for PYTHIA-like loops:
	    if (p->GetFirstDaughter() < p->GetId()) continue;
	    if (p->GetLastDaughter() - p->GetFirstDaughter() <=0)  continue;
	    //printf("analysis have found particle:");
	    //p->ls();
	    myParticles[nparticles]=i;
	    nparticles++;
	}
    } 

    if (nparticles<1) return;
*/
    
    HEPParticleList* stableChildren=new HEPParticleList;
    
    HEPParticleListIterator particles(plist);

    for (HEPParticle *part=particles.first(); part!=0; part=particles.next() ) {
	// exclude the ones which status is wrong.
	if (! part->Decays() ) {
	    MCDEBUGPART("Particle does not decay [status code = "<<part->GetStatus()<<"]",part)
	    continue;
	}
	// find all stable children;
	
	// at first - create the list of all daughters: intermediate and final
	HEPParticleList daughterlist;
	part->GetDaughterList(&daughterlist);
	if (daughterlist.empty()) {
	    
	    MCDEBUGPART("DAUGHTERLIST EMPTY",part);

	    //EVENT->ls();
	    printf("daughterlist is empty! decaying particle is:\n");
	    //part->ls();
	    exit(-1);
	}

	MCDEBUGPART("Particle OK ",part)
	
	// now dynamicaly append daughters' daughters etc.
	HEPParticleListIterator daughters(daughterlist);
	
	for (HEPParticle *d = daughters.first(); d!=0; d=daughters.next() ) {
	    
	    // only if status OK.
	    if (d->Decays() ) {
	    
		// but only if decay is not suppressed!
		if (Setup::IsSuppressed(d->GetPDGId())) continue;
		
		d->GetDaughterList(&daughterlist); // the list modifies itself in recurrece.
	    }
	}

	ndecayproducts=0;    
	stableChildren->clear();

	// select only the stable ones and fill decayproducts array:
	for (HEPParticle *d = daughters.first(); d!=0; d=daughters.next() ) {
	    if (d->IsStable() || (Setup::IsSuppressed(d->GetPDGId())) ){	

		stableChildren->push_back(d);
		//decayproducts[ndecayproducts]=d->GetId();
		ndecayproducts++;
		//d->ls();
	    }

	}

       // apply user's tree analysis:
	//stableChildren->ls();
	//printf("is empty:\n");
	//printf("%i\n",stableChildren->empty());
	//printf("stable children size:%i\n",stableChildren->size());
        if (userTreeAnalysis) {


	    HEPParticleList* x=stableChildren;
	    //printf("user tree analysis is %s\n",Setup::UserTreeAnalysis);
	    //printf("trying to apply user analysis. list is %i %i\n",stableChildren,x);
            //int userresult=(Setup::UserTreeAnalysis)(part,x);
            //int userresult=EnergyThresholdLimit(part,stableChildren);
	    Setup::UTA_particle=part;
	    Setup::UTA_partlist=stableChildren;
	    
	    userTreeAnalysis->Execute("Setup::UTA_particle,Setup::UTA_partlist,Setup::UTA_nparams, Setup::UTA_params");
	    
            if(stableChildren->empty()) continue;
	    //if(!stableChildren->size()==1) event->ls(); ***** need to look at extra entries.
        }

//	exit(-1);
    
//	TDecayMode *mode = TDecayMode::CheckMode( part->GetPDGId(),ndecayproducts,decayproducts,1, Setup::order_matters);
	TDecayMode *mode = TDecayMode::CheckMode( part->GetPDGId(),stableChildren,1, Setup::order_matters);
	// (Hint: no need to check if the routine returns valid pointer! )

	// Fill TDecayMode object with the data describing decay products:
	
	if (Setup::debug_mode) { 
	    if (mode->fill_histos)
		std::cout<<"###DM:";
	    else
		std::cout<<"###DM[NOHISTOS]:";
		
	    mode->ls();
	    std::cout<<std::endl;
	};
	
//	mode->Fill(ndecayproducts,decayproducts);
	mode->Fill(stableChildren,weight);

    }    
/*    
    // Now we process them one by one.
    
    // find all stable children:
    for (i=0; i<nparticles;i++) {
	nchildren=0;
	p=EVENT->GetParticle(myParticles[i]);

	// record all direct daughters:
	for (int j=p->GetFirstDaughter(); j<=p->GetLastDaughter();j++) {
	    children[nchildren]=j;
	    nchildren++;	    	
	}
	
	// and daughters' daughters...
	int k=0;
	while (k<nchildren) {
	    HEPParticle *d=EVENT->GetParticle(children[k]);

	    // if unstable - append to the end of list,
	    // so it will finally be processed...
	    if ( d->Decays() ) {
		// but only if decay is not suppressed!
		if (Setup::IsSuppressed(d->GetPDGId())) {
		    //printf("Decay Suppressed:");
		    //d->ls();
		} else {
		    for (int l=d->GetFirstDaughter();l<=d->GetLastDaughter();l++) {
			children[nchildren]=l;
			nchildren++;
		    }
		}
	    }
	    
	    k++;
	    
	} // end of while(k<nchildren)
	

	// We have a list of all decays products,
	// both: stable and unstable.
	//
	// we choose only final, stable products and fill
	// ndecayproducts table.
	
	ndecayproducts=0;
	for (int m=0; m<nchildren;m++){
	    HEPParticle *particle=EVENT->GetParticle(children[m]);
	    if (particle->IsStable() || (Setup::IsSuppressed(particle->GetPDGId())) ){	
		decayproducts[ndecayproducts]=particle->GetId();
		ndecayproducts++;
	    }
	}

	//printf("we look for decay mode with %i decay prods\n",ndecayproducts);

	// Let's find out whether we already have a TDecayMode object
	// that coresponds to our list of decay products.
	// If it did not exist, it will be created automatically for us.

	TDecayMode *mode = TDecayMode::CheckMode( p->GetPDGId(),ndecayproducts,decayproducts,1, Setup::order_matters);
	// (Hint: no need to check if the routine returns valid pointer! )

	// Fill TDecayMode object with the data describing decay products:
	
	mode->Fill(ndecayproducts,decayproducts);
    }
*/    

    if (Setup::user_event_analysis) {
	UserEventAnalysis &a=*(Setup::user_event_analysis);
	a.RestoreOriginalEvent(EVENT);
	Setup::EVENT=EVENT_ORIG;
    }


    delete stableChildren;

}

void PrintAnalysedEvent()
{
    if (Setup::user_event_analysis) {
	printf("event as modified by the user analysis code:\n");
	Setup::user_event_analysis->SavedEvent()->ls();
    }

}


void MC_FillUserHistogram(char *name, double value, double weight)
{
    TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
    assert(h);
    h->Fill(value,weight);

}