#include <CaloShape.h>
////////////////////////////////////////////////////////////
/// INITIALIZATIONS & GLOBAL PARAMETERS
////////////////////////////////////////////////////////////

CaloShape_parameters * CaloShape_parameters::_parameters = 0; 

/**
 * Set external parameters to default values
 */
void CaloShape_parameters::SetDefault(){
    return;
};
////////////////////////////////////////////////////////////
/// CLASS IMPLEMENTATION
////////////////////////////////////////////////////////////


void CaloShape::Init(){
    for(int ip=0; ip<23; ip++){
	qtot[ip]=0.;
	delta[ip]=0.;
	rms[ip]=0.;
	skew[ip]=0.;
	kurt[ip]=0.;
	qtrk[ip]=0.;
	mean[ip]=0.;
	nstrip[ip]=0;
    }
    qmax    = 0.;
    qtotmax = 0.;
    pmax    = -1;
}


bool CaloShape::Set(PamLevel2* l2, int ntr, int view, Bool_t usemechanicalalignement){
    if(!l2)return false;
    CaloTrkVar *c = l2->GetTrack(ntr)->GetCaloTrack();
    if(!c)return false;
    CaloLevel1* l1 = l2->GetCaloLevel1();
    if(!l1)return false;    
    Set(l1,c,view,usemechanicalalignement);
    return true;
};

bool CaloShape::Set(CaloLevel1* ca1, CaloTrkVar* ca2, int view, Bool_t usemechanicalalignement){

    Reset();
//	cout << " CaloEvent::Set()"<<endl; 

    if(!ca1)return false;
    CaloStrip cstrip;
    cstrip = CaloStrip(ca1,usemechanicalalignement);  

    for(int ih=0;ih<ca1->istrip;ih++){
	    
	int iv=-1;
	int ip=-1;
	int is=-1;
	ca1->DecodeEstrip(ih,iv,ip,is);

	// ----------------------
	// track intersection
	// ----------------------
	int is0 = -100;
	if(ca2) is0 = ca2->tibar[ip][iv] - 1;

	if(iv==view){
	    cstrip.Set(view,ip,is);
	    float fxy = 0.;
	    if( !view ) fxy = (float) (cstrip.GetX());
	    else        fxy = (float) (cstrip.GetY());
	    float fq = (float) (cstrip.GetE());		
	    nstrip[ip]++;
	    qtot[ip]  += fq;
	    delta[ip] += fq*fxy;
	    rms[ip]   += fq*fxy*fxy;
	    skew[ip]  += fq*fxy*fxy*fxy;
	    kurt[ip]  += fq*fxy*fxy*fxy*fxy;
	    if( is==is0 || is==is0-1 || is==is0+1 )qtrk[ip]+=fq;
	    if(fq>qmax)qmax = fq;
		
	}	
    }// end loop over hit strips

/* 	cout << "===============" <<endl; */
/* 	cout << qtot[4] << endl;  */
/* 	cout << delta[4] << endl;  */
/* 	cout << rms[4] << endl;  */
/* 	cout << skew[4] << endl;  */
/* 	cout << kurt[4] << endl;  */

    for(int ip=0; ip<22; ip++){

	float mu1t = 0.;
	float mu2t = 0.;
	float mu3t = 0.;
	float mu4t = 0.;
	if( qtot[ip]>0. ){

	    mu1t = delta[ip]/qtot[ip];
	    mu2t = rms[ip]/qtot[ip];
	    mu3t = skew[ip]/qtot[ip];
	    mu4t = kurt[ip]/qtot[ip];		

/* 		delta[ip] = mu1t; */
/* 		rms[ip] = mu2t - mu1t*mu1t; */
/* 		skew[ip] = mu3t - 3*mu1t*mu2t + 2*mu1t*mu1t*mu1t; */
/* 		kurt[ip] = mu4t - 4*mu3t*mu1t + 6*mu2t*mu1t*mu1t - 3*mu1t*mu1t*mu1t*mu1t; */

	    mean[ip] = mu1t;
	    // ----------------------
	    // track intersection
	    // ----------------------
	    float fxy0 = mu1t;
	    if(ca2)fxy0 = (float) (ca2->tbar[ip][view]);
	    //
	    delta[ip] = mu1t - fxy0;
	    rms[ip]   = mu2t - 2*mu1t*fxy0 + fxy0*fxy0;
	    skew[ip]  = mu3t - 3*mu2t*fxy0 + 3*mu1t*fxy0*fxy0 - fxy0*fxy0*fxy0;
	    kurt[ip]  = mu4t - 4*mu3t*fxy0 + 6*mu2t*fxy0*fxy0 - 4*mu1t*fxy0*fxy0*fxy0 + fxy0*fxy0*fxy0*fxy0;

	}

//	    if( nstrip[ip]>2 ){
	nstrip[22] += nstrip[ip];
	qtot[22]   += qtot[ip];
	delta[22]  += qtot[ip] * delta[ip];
	rms[22]    += qtot[ip] * rms[ip];
	skew[22]   += qtot[ip] * skew[ip];
	kurt[22]   += qtot[ip] * kurt[ip];
	mean[22]   += qtot[ip] * mean[ip];
	qtrk[22]   += qtrk[ip];
//	    }

	if(rms[ip]>0.)rms[ip] = sqrt( rms[ip] );
	else          rms[ip] = 0.;
	if(rms[ip]>0.)skew[ip] = skew[ip]/pow(rms[ip],3.);
//	if(rms[ip]>0.)skew[ip] = skew[ip];
	else          skew[ip] = 0.;
	if(rms[ip]>0.)kurt[ip] = kurt[ip]/pow(rms[ip],4.) - 3.;
//	if(rms[ip]>0.)kurt[ip] = kurt[ip]; 
	else          kurt[ip] = 0.;

	if( qtot[ip]>qtotmax ){
	    qtotmax = qtot[ip];
	    pmax = ip;
	}
	   
    }	

    if(qtot[22]>0.)mean[22] = mean[22]/qtot[22];
    else           mean[22] = 0.;
    if(qtot[22]>0.)delta[22] = delta[22]/qtot[22];
    else           delta[22] = 0.;
    if(rms[22]>0.)rms[22] = sqrt( rms[22]/qtot[22] );
    else          rms[22] = 0.;
    if(rms[22]>0.)skew[22] = skew[22]/( qtot[22]*pow(rms[22],3.) );
//    if(rms[22]>0.)skew[22] = skew[22]/( qtot[22] ) ;
    else          skew[22] = 0.;
    if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22]*pow(rms[22],4.) ) - 3.;
//    if(rms[22]>0.)kurt[22] = kurt[22]/( qtot[22] );
    else          kurt[22] = 0.;
	

/* 	cout << "===============" <<endl; */
/* 	cout << qtot[4]  << endl;  */
/* 	cout << delta[4] << endl;  */
/* 	cout << rms[4]   << endl;  */
/* 	cout << skew[4]  << endl;  */
/* 	cout << kurt[4]  << endl;  */

    return true;

};
///
void CaloShape::Print(){

    cout << endl;
//    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    cout << setw(3)<< "IP";
    cout << setw(5)<< "N";
    cout << setw(13)<< "Q";
    cout << setw(13)<< "DX";
    cout << setw(13)<< "RMS";
    cout << setw(13)<< "skew";
    cout << setw(13)<< "kurt";
    cout << endl;
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    for(int ip=0; ip<=22; ip++){
	cout << setw(3)<< ip;
	cout << setw(5)<< nstrip[ip];
	cout << setw(13)<< qtot[ip];
	cout << setw(13)<< delta[ip];
	cout << setw(13)<< rms[ip];
	cout << setw(13)<< skew[ip];
	cout << setw(13)<< kurt[ip];
	cout << endl;
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
    cout << endl;
}


ClassImp(CaloShape);
ClassImp(CaloShape_parameters);

