#include "TDecayMode.H"
#include "HEPEvent.H"
#include "TFile.h"
#include "TH1.h"
#include <iostream>

#include "Setup.H"

#include <stdlib.h>
#include <list>
#include <algorithm>

#include <math.h>

using namespace std;
ClassImp(TDecayMode)


int TDecayMode::NDecayModes=0;
long int TDecayMode::NFills=0;
TObjArray* TDecayMode::DecayModes=new TObjArray();
int MAX_MODES=200;

TDecayMode::TDecayMode():
mother_id(0),
num_of_daughters(0),
nentries(0),
sumw(0.0),
sumw2(0.0),
my_directory(0),
histograms(0),
fit_parameter(0)
{}

int TDecayMode::GetDaughter(int i)
{

    if (i<1 || i>num_of_daughters) {
	printf("ERROR: TDecayMode::GetDaughter(i=%i) out of range!\n",i);
	return 0; 
    }

    return daughters[i-1];
};

TDecayMode::TDecayMode(int motherid,HEPParticleList *daughterlist):
mother_id(motherid),
num_of_daughters(0),
nentries(0),
sumw(0.0),
sumw2(0.0),
my_directory(0),
fit_parameter(0)
{
//    printf("TDecayMode::TDecayMode\n");
//    HEPEvent &EVENT=*(Setup::EVENT);

    char buf[128];
    sprintf(buf,"%s => ",HEPParticle::GetParticleName(motherid));
        
    sprintf(latexname,"%s \\rightarrow ",HEPParticle::GetLaTeXName(motherid));

    //printf("latexname: %s\n",latexname);
/*
    // najpierw policzmy ile jest faktycznie stabilnych corek...
    num_of_daughters=0;
    for (int i=0; i<daughterlist_length;i++) {
	HEPParticle *p=EVENT.GetParticle(daughterlist[i]);
	if ( (p->IsStable()) || (Setup::IsSuppressed(p->GetPDGId() ) ) ) 
	    num_of_daughters++;    
    }
*/

    HEPParticleListIterator daughteritr(*daughterlist);
    for (HEPParticle *d=daughteritr.first(); d!=0; d=daughteritr.next()) {
	if ( (d->IsStable()) || (Setup::IsSuppressed(d->GetPDGId() ) ) ) 
	    num_of_daughters++;    
    }
    //printf("determined num_of_daughters=%i\n",num_of_daughters);
    // a teraz tworzymy liste kodow PDG samych stabilnych...
    int j=0;
    //printf("Num of real daughters:: %i\n",num_of_daughters);

    for (HEPParticle *d=daughteritr.first(); d!=0; d=daughteritr.next() ){
	if ( (d->IsStable()) || (Setup::IsSuppressed(d->GetPDGId()) ) ) {
	    daughters[j]=d->GetPDGId();
	    j++;
	    strcat(buf, d->GetParticleName());
	    strcat(buf," ");
	    strcat(latexname,d->GetLaTeXName());
	    strcat(latexname," ");
	    //printf("... %s\n",buf);
	}
    }

    //printf("  established name:%s\n",buf);
    SetName(buf);
    
    char varname[256];
    sprintf(varname,"DecayMode%03i",TDecayMode::NDecayModes);
    //printf("  established title:%s\n",varname);
    SetTitle(varname);
    gFile->cd();
//    my_directory=gFile->mkdir(varname,buf);
    my_directory=gFile->mkdir(buf,varname);
    //printf("Made directory: from:%s,%s\n",buf,varname);
    //my_directory->ls();
    
    my_directory->cd();
    my_directory->Append(this);
    //printf("#########ADDING MODE:%s#############\n",GetName());
    // now we'll create histograms:

    histograms=new TObjArray();

    if (num_of_daughters <= 7 && TDecayMode::NDecayModes < MAX_MODES) {
	fill_histos=1;
    } else {
	fill_histos=0;
    }

    if (fill_histos) {

    char hname[128];
    char htitle[128];
    // first: histograms of pairs:
    for (int i=0;i<num_of_daughters;i++) {
	for (int j=i+1;j<num_of_daughters;j++){
	
	   // histograms of mass:    
	   sprintf(hname,"hM1_%02i%02i",i,j);
	   sprintf(htitle,"Mass(1) of %s %s",HEPParticle::GetParticleName(daughters[i]),
	                                    HEPParticle::GetParticleName(daughters[j]));
	   //printf("adding histogram NAME:%s TITLE:%s\n",hname,htitle);
	   
	   int nbody=GetNumOfDaughters();
	   
	   TH1D *h=new TH1D(hname,htitle,
			    Setup::nbins  [nbody][2],
			    Setup::bin_min[nbody][2],
			    Setup::bin_max[nbody][2]);
	   
	   h->Sumw2();
	   //printf("Added histo (2):");
	   //h->ls();
	   histograms->Add(h);
	   //directory->Append(h);				    
	}
    }
    // now loop on a running list of existing ones and add further combination
    // up till all particles are there.
//printf("---BUILDING OTHER HISTOGRAMMES:\n");    


TH1D *h=0;
char hnumber[128];
char hId[10];
int kkk=0;
for (int hk=0;hk<=histograms->GetLast();hk++) {
    h=(TH1D*)histograms->At(hk);
    //printf("***[%i/%i]*** SCAN:",hk,histograms->GetLast());
    //h->ls();
    sscanf(h->GetName(),"%4s%s",hId,hnumber);
    int nextnum=(strlen(hnumber)/2);

    if (nextnum>num_of_daughters) {
    
    //printf("^^^NEXTNUM=%i > num_of_daughters=%i\n",nextnum,num_of_daughters);
    break;
    }

    // now we wanna get the last 2 digits and convert it to a number...
    int hptr=strlen(hnumber)-2;
    sscanf(&(hnumber[hptr]),"%i",&kkk);    
    //printf("This histogram has type [%s] and number[%s][%i] NEXT is %i+1 body, num_of_dayghters=%i\n",hId,hnumber,kkk,nextnum,num_of_daughters);
    if (nextnum>=num_of_daughters){
     //printf("Nothing more to add ...\n");
     break;
     }
     
     
    for (int i=nextnum;i<num_of_daughters;i++){ 

	   if (i<=kkk) {
	   //printf("KKK=%i is < than i=%i\n",kkk,i);
	   continue;
	   }
	   sprintf(hname,"%4s%s%02i",hId,hnumber,i);
	   sprintf(htitle,"%s %s",h->GetTitle(),HEPParticle::GetParticleName(daughters[i]));
	   //printf("NEW: NAME:%s TITLE:%s\n",hname,htitle);

	   int nbody=GetNumOfDaughters();


	   TH1D *h=new TH1D(hname,htitle,
				    Setup::nbins  [nbody][nextnum+1],
				    Setup::bin_min[nbody][nextnum+1],
				    Setup::bin_max[nbody][nextnum+1]);
	   h->Sumw2();
	   histograms->Add(h);
	   
    
    }

    } // endif (fill_histos);
//printf("============================\n");
}



};

TDecayMode::~TDecayMode()
{
  /**  for (int hk=0;hk<=histograms->GetLast();hk++) {
    delete histograms->At(hk);
  }**/
  if (histograms) delete histograms;
};


int ComparePDGs(HEPParticle* p1, HEPParticle* p2)
{
    double pdg1=(double)p1->GetPDGId();
    if (pdg1<0) pdg1=-pdg1+0.5;
    double pdg2=(double)p2->GetPDGId();
    if (pdg2<0) pdg2=-pdg2+0.5;

    //printf("Comparing particles: pdg1=%f pdg2=%f (id1=%i id2=%i)\n",pdg1,pdg2,p1->GetId(),p2->GetId());

    if (pdg1>pdg2) return 1;
    //else if (pdg1<pdg2) return -1;
    else return 0;

}


int CompareParticlePDGs(const void* a1, const void* a2)
{
    HEPEvent &EVENT=*(Setup::EVENT);


    int *idx1=(int*)(a1);
    int *idx2=(int*)(a2);
    
    HEPParticle *p1=EVENT.GetParticle(*idx1);
    HEPParticle *p2=EVENT.GetParticle(*idx2);

    return ComparePDGs(p1,p2);
}

TDecayMode* TDecayMode::CheckMode(int motherid,HEPParticleList *daughterlist, int autoadd, int order_matters)
{
    //printf("CheckMode:HepParticleList ndecayproducts=%i\n",daughterlist->size());
    //HEPParticleListIterator itr(*daughterlist);
    //for (HEPParticle *part=itr.first(); part!=0; part=itr.next() ) part->ls();
    //printf("\n---\n");

//    HEPEvent &EVENT=*(Setup::EVENT);
    
    if (!order_matters) {
	std::list<HEPParticle*> *l = (std::list<HEPParticle*>*) daughterlist->GetList();
	
	l->sort(ComparePDGs);
//	sort((*l).begin(),(*l).end(),ComparePDGs);

//    // we must sort decayproducts using qsort
//    qsort(decayproducts,ndecayproducts,sizeof(int),CompareParticlePDGs);
    }

    //printf("AFTER SORT : CheckMode:HepParticleList ndecayproducts=%i\n",daughterlist->size());
    //HEPParticleListIterator itr2(*daughterlist);
    //for (HEPParticle *part=itr2.first(); part!=0; part=itr2.next() ) part->ls();
    //printf("\n---\n");


    char modename[128];
    
    //printf("TDecayMode::CheckMode(%i, %i)\n",ndecayproducts,decayproducts);


//    if (!ndecayproducts || ! decayproducts) {    
    if (daughterlist->empty()) {    
	printf("ERROR in TDecayMode::CheckMode. No decay products!\n");
	exit(-1);
    }

    sprintf(modename,"%s => ",HEPParticle::GetParticleName(motherid));

    HEPParticleListIterator daughteritr(*daughterlist);


    for (HEPParticle *d=daughteritr.first(); d!=0; d=daughteritr.next() ) {
    
	strcat(modename,d->GetParticleName());
	strcat(modename," ");
    }
    
    //printf("MODE NAME:%s\n",modename);
    
    //printf("looking for object:[%s]\n",modename);
    TDecayMode *dm=(TDecayMode*)(TDecayMode::DecayModes)->FindObject(modename);    

    //printf("TDecayMode tried to find:%i\n",dm);
    //printf("Found: %i \n",dm);
    if ((!dm) && autoadd ) {
	//printf("needs to add mode...\n");
	//dm= new TDecayMode(motherid,ndecayproducts,decayproducts);
	dm= new TDecayMode(motherid,daughterlist);
	//printf("Created Decay mode :%s\n",dm->GetName());
	TDecayMode::DecayModes->AddLast(dm);
	TDecayMode::NDecayModes++;
	//printf("added mode:");
	//dm->ls();
    } else {
	//printf("found mode:");
	//dm->ls();
    }
    //printf("finally: MODE=%i\n",dm);
    //printf("**********************\n");
    
    return dm;
};


void TDecayMode::Fill(HEPParticleList *daughterlist, double weight)
{
    //printf("TDecayMode::Fill():");

    //printf("FILL: ndecayproducts=%i\n",daughterlist->size());
    //HEPParticleListIterator itr(*daughterlist);
    //for (HEPParticle *part=itr.first(); part!=0; part=itr.next() ) part->ls();
    //printf("\n---\n");

    
//    HEPEvent &EVENT=*(Setup::EVENT);

    nentries++;
    sumw+=weight;
    sumw2+=(weight*weight);
    TDecayMode::NFills++;

    if (!fill_histos) return;

    int nhist=histograms->GetLast();

    if (nhist<0){    
     printf ("WARNING: TDecayMode::Fill for mode %s -> nhist=%i\n",GetName(),nhist);
     histograms->ls();
     printf("********\n");
    }

    // rewrite list to array:
    HEPParticle* decayproducts[100];
    HEPParticleListIterator daughteritr(*daughterlist);
    int ind=0;
    for (HEPParticle *d=daughteritr.first(); d!=0;d=daughteritr.next()) {
	decayproducts[ind]=d;
	ind++;
    }


    for (int i=0; i<=nhist;i++) {

    TH1D *h=(TH1D*)histograms->At(i);
    //printf("FILLING:");
    //h->ls();
    
    double E=0;
    double PX=0;
    double PY=0;
    double PZ=0;
    
    char mod;
    int pwr;
    char num[128];
    //                     h M 2 _ 010203
    sscanf(h->GetName(),"h%c%1i_%s",&mod,&pwr,num);
    
    //printf("This is histo of %c in power %i for daughters %s\n",mod,pwr,num);


    if (mod=='M') {
	int nbody=strlen(num)/2;
	for (int j=0;j<nbody;j++) {
	    int idx=0;
	    sscanf(&(num[2*j]),"%2i",&idx);
	    HEPParticle *p=decayproducts[idx];
	//printf("inquiring particle:");
	    //p->ls();
	    E+=p->GetE();
	    PX+=p->GetPx();
	    PY+=p->GetPy();
	    PZ+=p->GetPz();
	}
	double P2=PX*PX+PY*PY+PZ*PZ;
	double M=sqrt(E*E-P2);
	if (pwr>1) {
	    double m1=M;
	    for(int mm=1; mm<pwr;mm++) M=M*m1;
	}
	h->Fill(M,weight);
	    
    } else {
	printf("HISTOGRAM WITH UNKNOWN MODE : %s\n",h->GetName());
    }


   }
}

void TDecayMode::ls(char *option)
{
printf("TDecayMode:[%s] :%li entries\n",GetName(),GetNEntries());
}


void TDecayMode::SetMyDirectory(TDirectory *d)
{
//    cout<<"SetMyDirectory called:"<<d<<"\n";
    this->my_directory=d;
}

TDirectory* TDecayMode::GetMyDirectory()
{
//    cout<<"GetMyDirectory:"<<this->my_directory<<"\n";
    return this->my_directory;
}