#include <math.h>
#include <new.h>
#include <iostream.h>
#include "archparam.h"
#include "pisteetoile.h"
 
TransFuncElec PisteEtoile::Trapani;
FormePulse    PisteEtoile::FPulse(Trapani);

// D. Yvon, CE Saclay, DAPNIA/SPP, 08/99

PisteEtoile::PisteEtoile(short NoPhotDiod) {
	PasTempsEch=archParam.acq.perEch;			
		//BUGG doit changer	!!!!! Doit etre extrait d'un serveur de param d'acquisition
	NoPhDiod=NoPhotDiod;
	DiodVivante=true;
	
	SSTnbPasFit= (int) (2*TFExcursion/TPasTFit);
	
	AmplNorm=0.;
	for(int i=0; i<SSTLongIndexEvt()*10; i++){
		AmplNorm+=FPulse.PulseShape(i*PasTempsEch/10.)*FPulse.PulseShape(i*PasTempsEch/10.);
	}
	AmplNorm/=10.;
	indexDebutPiste=0;
	
	PhDiodArray=new(nothrow) int[PhDiodTabLong];
	PhDiodPiste=new(nothrow) double[PhDiodTabLong];
	Ampl=new(nothrow) double[SSTnbPasFit];
	X2=new(nothrow) double[SSTnbPasFit];
	
	if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
		cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
		exit(-1);
	}
}

PisteEtoile& PisteEtoile::operator =(PisteEtoile const & x) {
	PasTempsEch=x.PasTempsEch;		
	NoPhDiod=x.NoPhDiod;
	DiodVivante=x.DiodVivante;
	indexDebutPiste=x.indexDebutPiste;
	DiodVivante=x.DiodVivante;
	AmplNorm=x.AmplNorm;

	
	for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
	for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
	
	LastEtoile=x.LastEtoile;
	
	 gachette=x.gachette;
	 Vrms=x.Vrms;
	 MaxIndex=x.MaxIndex;			// Index de l'amplitude Max en analyse Trig()
	 AmplMax=x.AmplMax;
	 
	 AmplRes=x.AmplRes;
	 TIndex=x.TIndex;
	 X2Min=x.X2Min;
	 SuccesAFine=x.SuccesAFine;
	 AmplNorm=x.AmplNorm;
	 
	 return *this;
}

PisteEtoile::PisteEtoile (PisteEtoile const & x) {
	PasTempsEch=x.PasTempsEch;		
	NoPhDiod=x.NoPhDiod;
	DiodVivante=x.DiodVivante;
	indexDebutPiste=x.indexDebutPiste;
	DiodVivante=x.DiodVivante;
	AmplNorm=x.AmplNorm;
		
	PhDiodArray=new(nothrow) int[PhDiodTabLong];
	PhDiodPiste=new(nothrow) double[PhDiodTabLong];
	Ampl=new(nothrow) double[SSTnbPasFit];
	X2=new(nothrow) double[SSTnbPasFit];
	
	
	if((X2==0)||(Ampl==0)||(PhDiodPiste==0)||(PhDiodArray==0)) {
		cerr<<"Erreur a la reservation de memoire dans constructeur pisteetoile"<<endl;
		exit(-1);
	}
	
	for(int j=0; j<PhDiodTabLong;j++) PhDiodArray[j]=x.PhDiodArray[j];
	for(int j=0; j<PhDiodTabLong;j++) PhDiodPiste[j]=x.PhDiodPiste[j];
	
	LastEtoile=x.LastEtoile;
	
	 gachette=x.gachette;
	 Vrms=x.Vrms;
	 MaxIndex=x.MaxIndex;			// Index de l'amplitude Max en analyse Trig()
	 AmplMax=x.AmplMax;
	 
	 AmplRes=x.AmplRes;
	 TIndex=x.TIndex;
	 X2Min=x.X2Min;
	 SuccesAFine=x.SuccesAFine;
	 AmplNorm=x.AmplNorm;

}



PisteEtoile::~PisteEtoile() {
	delete[] X2;
	delete[] Ampl;
	delete[] PhDiodPiste;
	delete[] PhDiodArray;
	
}

void PisteEtoile::fill(int* pData, int InDebutPiste, int nData) {
//	if(!isAlive()) return;
	indexDebutPiste=InDebutPiste - PhDiodTabLong + nData;
	int* PhDiod=PhDiodArray + PhDiodTabLong - nData;
	int* pDataC=pData;
	for(int k=0;k<nData;k++,pDataC++, PhDiod++) (*PhDiod)=(*pDataC);
	return;
}

void PisteEtoile::push(int* pData, int nbdata) {
	int* PPhDiod=PhDiodArray;
	int* PPhDiod2=PhDiodArray+nbdata;
	int* pDataC=pData;
	
	for(int i=0; i<PhDiodTabLong-nbdata; i++, PPhDiod++, PPhDiod2++) *PPhDiod=(*PPhDiod2);
	for(int i=0; i<nbdata; i++, PPhDiod++,pDataC++) *PPhDiod=(*pDataC);
	indexDebutPiste+=nbdata; 
	return;
}

bool PisteEtoile::traque() {
	if(!isAlive()) return false;		//La diode est-elle morte?
	
	bool testEtoile=trig();
	if(testEtoile) {			
		//On a declenché sur une Etoile
		if(archParam.sst.analFine) {
			demasque();
// Un second Filtre est possible
			// if(!(fabs(AmplRes)>SeuilGachette*Vrms)&&SuccesAFine) testEtoile=false;
			// Test rajouté pour filtrer le bruit en escalier des voies SST
			// Le fit de forme tasse le bruit et permet un meilleur filtrage.
		}

// Troisième filtre possible:
		// On tue les evenements trop larges
		if (widthHalfMax>(FPulse.FWHM()/PasTempsEch+2.)) testEtoile=false;
	
		if(testEtoile) {
		// On emballe une etoile
			LastEtoile.InpCurrent= AmplRes/FPulse.GainElectrique();
			LastEtoile.InpCurFin= AmplFin/FPulse.GainElectrique();
			LastEtoile.TEchan=TIndex;
			LastEtoile.TFin=Tfin;
			LastEtoile.InpCurrNoisRMS=Vrms/FPulse.GainElectrique();
			LastEtoile.X2Calc=X2Min;
			LastEtoile.NoDiode=NoPhDiod;
			LastEtoile.FitForme=SuccesAFine;
			LastEtoile.LargMiHauteur=widthHalfMax;
			
/*			cout<<setprecision(9);
			cout<<TIndex<<" "<<Tfin<<" "<<AmplRes<<" "<<AmplFin<<endl;
			cout<<TIndex-Tfin<<" "<<(AmplRes-AmplFin)/AmplFin<<endl;
*/
		}
	}
	// On expedie	
	return testEtoile;
}

bool PisteEtoile::trig() {	
	MaxIndex=-1;
	AmplMax=0.;
	gachette=false;

// Copie les données brutes sur la piste. Les met en Volt ADC 12 Bits-10 Volts
	for(int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]=PhDiodArray[i]*SSTLsbVal;

// Calcule Valeur Moy du prepulse
	double Moy=0.;
	for(int i=0; i<Prepulselong; i++) Moy+=PhDiodPiste[i];
	Moy/=Prepulselong;

// On soustrait la valeur Moy de la piste
	for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=Moy;
	
	double pente=0.;
	if (archParam.sst.soustPente==true) {
	//	On calcule la pente du prepulse
		double milprepulse=(Prepulselong-1.)/2.;
		double SomCarre=0.;
		for(int i=0; i<Prepulselong; i++) {
			SomCarre+= (i-milprepulse)*(i-milprepulse);
			pente+= (i-milprepulse)*PhDiodPiste[i];
		}
		pente=pente/SomCarre;

	// On soustrait la pente sur la piste
		for (int i=0; i<PhDiodTabLong; i++) PhDiodPiste[i]-=(i-milprepulse)*pente;
	}
	
// On calcule le Vrms sur le prepulse
	Vrms=0.;
	for (int i=0; i<Prepulselong; i++) {
		Vrms+=PhDiodPiste[i]*PhDiodPiste[i];
	}
	Vrms=sqrt(Vrms/Prepulselong);
	if(Vrms<SSTLsbVal) Vrms=SSTLsbVal;		// Pour eviter les VRMS=0.
	
// Declenche-t-on?
	gachette=false;
	short iSeuil=0;
	double* pPiste =PhDiodPiste+Prepulselong-1;
	
// On demande deux échantillons à seuilgachette Vrms
	for (int i=0; i<(PhDiodTabLong-Prepulselong); i++, pPiste++) {
		if (fabs(*pPiste)>(SeuilGachette*Vrms)) 
			if (fabs(*(pPiste+1))>SeuilGachette*Vrms) {
				gachette=true;
				iSeuil=Prepulselong-1+i;
				break;
			}		
	}

	if (gachette==true) {
// On trouve le maximum
		MaxIndex=iSeuil;
		AmplMax=fabs(PhDiodPiste[iSeuil]);
		for (int i=iSeuil; i<PhDiodTabLong; i++) {
			if (fabs(PhDiodPiste[i])>AmplMax) {
				MaxIndex=i;
				AmplMax=fabs(PhDiodPiste[i]);
			}
		}
		if ((MaxIndex>=PhDiodTabLong-4)||(MaxIndex<Prepulselong+4)) gachette=false;
// Cela sera traite a l'iteration suivante
// Ou ca a du etre traite a l'iteration precedente	
	}
	if (gachette==true) {
// On calcule la largeur a mi hauteur
		int HalfMaxInf;
		for(HalfMaxInf=MaxIndex;HalfMaxInf>Prepulselong;HalfMaxInf--)
			if(fabs(PhDiodPiste[HalfMaxInf])<(AmplMax/2.)) break;
				int HalfMaxSup;	
		for(HalfMaxSup=MaxIndex;HalfMaxSup<PhDiodTabLong;HalfMaxSup++)
			if(fabs(PhDiodPiste[HalfMaxSup])<(AmplMax/2.)) break;
		widthHalfMax=HalfMaxSup-HalfMaxInf;
		
/*
// On integre sur la larg a mi hauteur
		integ=0.;
		if((HalfMaxSup<PhDiodTabLong)&&(HalfMaxInf>Prepulselong))
			for(int i=HalfMaxInf; i<=HalfMaxSup; i++) integ+=PhDiodPiste[i];
*/


// On fitte une parabole autour du maximum
		double tmaxfit;
		if (widthHalfMax>1) {
			double ymax=PhDiodPiste[MaxIndex];
			double yplus=PhDiodPiste[MaxIndex+1];
			double ymoins=PhDiodPiste[MaxIndex-1];
	// Estimation du temps
			if(!(ymax==yplus)) {
				double Y=(yplus-ymax)/(ymax-ymoins);
				tmaxfit=0.5*(1.+Y)/(1-Y);
			}
			else tmaxfit=0.5;
			
	// Estimation de l'amplitude
			double dtplus2=(1.+tmaxfit)*(1.+tmaxfit);
			double dtmoins2=(1.-tmaxfit)*(1.-tmaxfit);
			if(dtplus2==dtmoins2) AmplRes=PhDiodPiste[MaxIndex];
			else AmplRes= (yplus*dtplus2-ymoins*dtmoins2) / (dtplus2-dtmoins2);
		}
		
//	On stocke les resultats
		TIndex=indexDebutPiste+MaxIndex+tmaxfit;
		X2Min=0.;
		SuccesAFine=false;
	}
	return gachette;	
}

void PisteEtoile::demasque() {
// Maintenant que l'on se doute que il y a une etoile 
// On veut connaitre son flux et son temps de passage
// Nécessite de connaitre la vitesse de rotation du ballon.

		// BUGG devra etre extrait des headers acquisition. Ca peut changer!!!!!!
			
	int IndexDebutEvt=MaxIndex-SSTLargEvtIndex();
	
	pFitArray0= PhDiodPiste+IndexDebutEvt; 
		// Pointeur origine du tableau d'echantillon sur lequel sera calculé le fit
	double DeltaTemps=0.;
	
// Pour chaque decalage en temps, on estime une amplitude et calcule le X2
	for(int k=0; k<SSTnbPasFit; k++) {
		DeltaTemps=TFExcursion-k*TPasTFit;
		pFitArray=pFitArray0;
		double Som=0.;
		for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++)	{
//			cout<<(*pFitArray)<<'\t'<<FPulse.PulseShape(i*PasTempsEch-DeltaTemps)<<endl;
			Som=Som+FPulse.PulseShape(i*PasTempsEch-DeltaTemps)*(*pFitArray);
		}
		Ampl[k]=Som/AmplNorm;
		X2[k]=0.;
		double X2Som=0.; double dum=0.;pFitArray=pFitArray0;
		for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) {
			dum=(*pFitArray)-Ampl[k]*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); 	// Différence
			X2Som+=dum*dum;															// Au carré
		}				
		X2[k]=X2Som/(SSTLongIndexEvt()*Vrms*Vrms);
	}
	
// On recherche le minimum du X2
	X2Min=1.e99;
	int X2MinIndex=-1;
	for(int k=0; k<SSTnbPasFit; k++) { 
		if (X2[k]<X2Min) {
			X2Min=X2[k];
			X2MinIndex=k;
		}
	}
	
	if((X2MinIndex>0)&&(X2MinIndex<(SSTnbPasFit-1))) {
	
// On fitte une parabole sur les trois X2 qui encadrent le min
		double Di;	// Index (non entier) du minimum du X2
		double X=(X2[X2MinIndex+1]-X2Min)/(X2Min-X2[X2MinIndex-1]);		
		Di=0.5*(X+1)/(X-1);

// Le X2 est mimimum pour le n° de pas de fit X2MinIndex+Di 
// On calcule l'amplitude estimee a ce temps
		double Som1=0.; pFitArray=pFitArray0;
		DeltaTemps=TFExcursion-(X2MinIndex-Di)*TPasTFit;
		for(int i=0; i<SSTLongIndexEvt(); i++, pFitArray++) 
				Som1+=(*pFitArray)*FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
		AmplFin=Som1/AmplNorm;
			


#ifdef PEtoileDebug
		double *PulseTrak=new double [SSTLongIndexEvt()];
#endif
// On calcule le X2 correspondant
		double X2Som1=0.; double dummy=0.;pFitArray=pFitArray0;
		for(int i=0;i<SSTLongIndexEvt();i++,pFitArray++) {	
		
#ifdef PEtoileDebug
			PulseTrak[i]=FPulse.PulseShape(i*PasTempsEch-DeltaTemps);
#endif
			dummy=(*pFitArray)-AmplFin*FPulse.PulseShape(i*PasTempsEch-DeltaTemps); 	// Différence
			X2Som1+=dummy*dummy;														// Au carré
		}
						
		X2Min=X2Som1/(SSTLongIndexEvt()*Vrms*Vrms);
		
		Tfin= indexDebutPiste+IndexDebutEvt+DeltaTemps/PasTempsEch;
			//Ouff!!!!
/*			
		cout<<setprecision(9)<<endl;			
		cout<<Tfin<<" "<<indexDebutPiste<<" "<<IndexDebutEvt<<" "<<DeltaTemps/PasTempsEch<<endl;
*/		
		SuccesAFine=true;
		
#ifdef PEtoileDebug
		delete[] PulseTrak;
#endif
	}
	else {				// Echec, c'est tres bizarre! On doit etre dans le bruit
		X2Min=0.;		
	}
	return;
} 

int  PisteEtoile::SSTLargEvtIndex() {  		
// Demi-largeur a 3 sigma de l'impulsion 
	int dum=(int) (15.e-3/PasTempsEch+1.); 
	return dum;
}

int  PisteEtoile::SSTLongIndexEvt() {		
// Nombre d'echantillon sur lequel est calculé le fit de forme
	int temp=2*SSTLargEvtIndex()+1;
	return temp;
}