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

CaloHough_parameters * CaloHough_parameters::_parameters = 0; 

/**
 * Set external parameters to default values
 */
void CaloHough_parameters::SetDefault(){

    cout << "CaloHough --> SET DEFAULT PARAMETERS" << endl;

    isimu = 0;

//     Q_UP = 4.;// mip
//     Q_DN = 0.;
//     M_UP = 3;
//     M_DN = 0;
    Q_UP = 10000.;
    Q_DN = 100.;
    M_UP = 3000;
    M_DN = 0;

    h_qtot[0] = new TH2F("houghqtot_0","Hit pattern",96,-0.5,95.5,22,-0.5,21.5);
    h_qtot[1] = new TH2F("houghqtot_1","Hit pattern",96,-0.5,95.5,22,-0.5,21.5);

    int n1 = 150;
    int n2 = 150;

    h_par[0]  = new TH2F("houghpar_0","par1 vs par2",n1,-40.,40.,n2,-10.,10.);
    h_par[1]  = new TH2F("houghpar_1","par1 vs par2",n1,-40.,40.,n2,-10.,10.);

    h_par_fine[0]  = new TH2F("houghparfine_0","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.);
    h_par_fine[1]  = new TH2F("houghparfine_1","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.);

    h_par1[0]  = new TH1F("houghpar1_0","par1",n1,-10.,10.);
    h_par1[1]  = new TH1F("houghpar1_1","par1",n1,-10.,10.);

    h_par2[0]  = new TH1F("houghpar2_0","par2",n2,-40.,40.);
    h_par2[1]  = new TH1F("houghpar2_1","par2",n2,-40.,40.);
};
////////////////////////////////////////////////////////////
/// CLASS IMPLEMENTATION
////////////////////////////////////////////////////////////
/**
 * Reset event
 * it has to be called for each event
 */
void CaloHough::Reset(){
    return;
}

/**
 * \brief Set event 
 */
bool CaloHough::Set(CaloLevel1 *cl1){

//    bool debug   = CaloHough_parameters::Get()->debug;    
    float Q_UP    = CaloHough_parameters::Get()->Q_UP;
    int   M_UP    = CaloHough_parameters::Get()->M_UP;
    float Q_DN    = CaloHough_parameters::Get()->Q_DN;
    int   M_DN    = CaloHough_parameters::Get()->M_DN;

    cout << "bool CaloHough::Set(CaloLevel1*)"<<endl;

    if(!cl1)return false;

    TH2F* h_qtot[2];
    TH2F* h_par[2];
    TH2F* h_par_fine[2];
    TH1F* h_par1[2];
    TH1F* h_par2[2];
    for(int iv=0; iv<2; iv++){
	h_qtot[iv] = CaloHough_parameters::Get()->h_qtot[iv];
	h_par[iv]  = CaloHough_parameters::Get()->h_par[iv];
	h_par_fine[iv]  = CaloHough_parameters::Get()->h_par_fine[iv];
	h_par1[iv] = CaloHough_parameters::Get()->h_par1[iv];
	h_par2[iv] = CaloHough_parameters::Get()->h_par2[iv];
	if(!h_qtot[iv] || !h_par[iv] || !h_par1[iv]  || !h_par2[iv] || !h_par_fine[iv])return false;
	h_qtot[iv]->Reset();
	h_par[iv]->Reset();
	h_par_fine[iv]->Reset();
	h_par1[iv]->Reset();
	h_par2[iv]->Reset();
    }

    // -----------------------------------
    // set charge collected in each plane
    // -----------------------------------
    //
    CaloStrip cstrip;
    int isimu = CaloHough_parameters::Get()->isimu;
    if(isimu==0)cstrip = CaloStrip(cl1,false);
    if(isimu==1)cstrip = CaloStrip(cl1,true);
    //
    vector<int>   ids[22][2];  //id strip 0-95
    vector<float>  xy[22][2];  //cluster coordinate
    int islast   = -1;
    int iplast   = -1;
    int mult     = 0;
    float qsum   = 0.;
    float xyqsum = 0.;
    float issum  = 0.;
    int ncluster[] = {0,0};
    for(int ih=0;ih<cl1->istrip;ih++){
	//
	int iv=-1;
	int ip=-1;
	int is=-1;
	cl1->DecodeEstrip(ih,iv,ip,is);
	float xyz[3];
	cstrip.Set(iv,ip,is);
//	cout << " * "<<setw(3)<<iv<<setw(3)<<ip<<setw(3)<<is<<endl;
	xyz[0] = cstrip.GetX();
	xyz[1] = cstrip.GetY();
	xyz[2] = cstrip.GetZ();
	//
//      ...................................... metodo 2
	if( is != islast+1 || ip != iplast ){ // new cluster
	    if( 
		qsum <  Q_UP && 
		qsum >  Q_DN && 
		mult <  M_UP && 
		mult >  M_DN && 
		true){         // store the previous cluster
		xy[ip][iv].push_back( xyqsum/qsum );
		ids[ip][iv].push_back( (int)(0.5+issum/qsum) );		
		h_qtot[iv]->Fill( issum/qsum, (float)(21-ip), qsum);		
		ncluster[iv]++;
		cout << " >> "<<qsum<<endl;
	    }// ... and reset counters	    
	    mult = 0;
	    qsum = 0.;
	    xyqsum = 0.;
	    issum  = 0.;  
	}

	mult++;
	qsum   +=             cstrip.GetE();
	xyqsum += xyz[iv]   * cstrip.GetE();
	issum  += (float)is * cstrip.GetE();
	islast  = is;
	iplast  = ip;	    

//      ...................................... metodo 1
// 	if( cstrip.GetE() < Q_UP ){
// 	    h_qtot[iv]->Fill( (float)is, (float)(21-ip), cstrip.GetE());
// 	    ids[ip][iv].push_back(is);
// 	}

    }
    cout << "n.clusters X "<<ncluster[0]<<endl;
    cout << "n.clusters Y "<<ncluster[0]<<endl;
    //
//    cstrip.Set(0,0,0);
//    cstrip.Set(0,21,0);
    cstrip.Set(0,11,0);
    float zref = cstrip.GetZ();//reference plane
    //
    for(int iv=0; iv<2; iv++){
	for(int ip1=0; ip1<22-1; ip1++){
//	    for(int ip2=ip1+1; ip2<22; ip2++){
	    for(int ip2=ip1+1; ip2<min(22,ip1+3); ip2++){
//		cout << "ip1 "<<ip1<<" "<<ids[ip1][iv].size()<<endl;
//		cout << "ip2 "<<ip2<<" "<<ids[ip2][iv].size()<<endl;
		for(int is1=0; is1<ids[ip1][iv].size(); is1++){
		    for(int is2=0; is2<ids[ip2][iv].size(); is2++){

//      ...................................... metodo 1
// 			int ids1 = ids[ip1][iv][is1];
// 			int ids2 = ids[ip2][iv][is2];
// 			float xyz1[3];
// 			float xyz2[3];
// 			cstrip.Set(iv,ip1,ids1);
// 			xyz1[0] = cstrip.GetX();
// 			xyz1[1] = cstrip.GetY();
// 			xyz1[2] = cstrip.GetZ();
// 			cstrip.Set(iv,ip2,ids2);
// 			xyz2[0] = cstrip.GetX();
// 			xyz2[1] = cstrip.GetY();
// 			xyz2[2] = cstrip.GetZ();

//      ...................................... metodo 2
 			int ids1 = ids[ip1][iv][is1];
 			int ids2 = ids[ip2][iv][is2];
 			float xyz1[3];
 			float xyz2[3];
 			cstrip.Set(iv,ip1,is1);
			xyz1[iv] = xy[ip1][iv][is1];
			xyz1[2]  = cstrip.GetZ();
			cstrip.Set(iv,ip2,ids2);
			xyz2[iv] = xy[ip2][iv][is2];
			xyz2[2]  = cstrip.GetZ();
			//
			float par1 = 0.;
			if(fabs(xyz1[2]-xyz2[2])>0.)par1 = (xyz1[iv]-xyz2[iv])/(xyz1[2]-xyz2[2]);
			float par2 = xyz1[iv] - par1*xyz1[2];
			par2 = par1*zref + par2;
			h_par[iv]->Fill(par2,par1); // parameter space
			h_par_fine[iv]->Fill(par2,par1); // parameter space
			h_par1[iv]->Fill(par1);
			h_par2[iv]->Fill(par2);
//			cout << setw(15)<<par1<< setw(15)<<par2<<endl;			
		    }
		}
	    }
	}
    //
	for( int ibx=1; ibx <= h_par[iv]->GetXaxis()->GetNbins(); ibx++){
	    for( int iby=1; iby <= h_par[iv]->GetYaxis()->GetNbins(); iby++){
		// get bin content
		float content = h_par[iv]->GetBinContent(ibx,iby);
		if( content > 6){
 		    double xmi = h_par[iv]->GetXaxis()->GetBinLowEdge(ibx);
 		    double xma = h_par[iv]->GetXaxis()->GetBinUpEdge(ibx);
 		    double ymi = h_par[iv]->GetYaxis()->GetBinLowEdge(iby);
 		    double yma = h_par[iv]->GetYaxis()->GetBinUpEdge(iby);
		    h_par_fine[iv]->GetXaxis()->SetRangeUser(xmi,xma);
		    h_par_fine[iv]->GetYaxis()->SetRangeUser(ymi,yma);
		    double par1 = h_par_fine[iv]->GetMean(1);
		    double par2 = h_par_fine[iv]->GetMean(2);
		    cout << iv << " <><><><><><><>< "<<setw(15)<<par1<<setw(15)<<par2<<endl;
		}	    
	    }
	}
	
	//
	for(int ip=0; ip<22; ip++){
	    ids[ip][iv].clear();
	    xy[ip][iv].clear();
	}	
    }//end loop on views

    return true;
}
ClassImp(CaloHough);
ClassImp(CaloHough_parameters);