// the function computes distances from many points to the parabola described by the input values a,b,c (using some global variables and a numetrical minimisation algorithm minDist).



double PRF_CHI_ev (const double *x){
    //cout <<"############## FUNCTION CALL ##########"<<endl;
    cnt++;
    double a = x[0];
    double b = x[1];
    double c = x[2];
    double chi = 0;
    int counter = 0;
    double error;
    bool sign;
    double sig,sigdn,sigup;
    double Localqf[70][10];
    double x_closest,y_closest,x_pad,y_pad;
    double distance, exp_distance;
    
    TGraphErrors * SingleEvPRF = new TGraphErrors();
    distandqfrac PadPair[70][10];
    //cout<<"---------------------------------------------------------"<<endl;
    //cout<<" a: "<<a<<" b: "<<b<<" c: "<<c<<endl;
    if (x[1]!=x[1] || x[0]!=x[0] || x[2]!=x[2]){
        cout<< "############## FUNCTION CALL ERROR ##########"<<endl;
        return 1.e13;
    }
    //cout<<func->GetParameter(0)<<" "<<func->GetParameter(1)<<" "<<func->GetParameter(2)<<" "<<func->GetParameter(3)<<endl;
    for(int j=0;j<clusters;j++){
        for(int k=0;k<multiplicity[j];k++){

            
            //cout<<chargemap->Interpolate(x_pos[padx[j][k]],y_pos[pady[j][k]])<<endl;
            //cout<<j<<" "<<k<<" "<<x_pos[padx[j][k]]<<" "<<y_pos[pady[j][k]]<<" "<<charge[j][k]<<"\t"<<a<<" "<<b<<" "<<c<<endl;
            thisSeg->SetParabola(x[0],x[1],x[2]);
            para->SetParameter(2,x[0]);
            para->SetParameter(1,x[1]);
            para->SetParameter(0,x[2]);
            PadPair[j][k] = ChargeFraction_real(x_pos[padx[j][k]],y_pos[pady[j][k]],charge[j][k],thisSeg);
            Localqf[j][k] = PadPair[j][k].qfrac;
            

            x_pad = x_pos[padx[j][k]];
            y_pad = y_pos[pady[j][k]];

            x_closest = minDist(a,b,c,x_pad,y_pad);
            y_closest = x_closest*x_closest*a + x_closest*b + c;
            distance = sqrt( (x_closest-x_pad)*(x_closest-x_pad) + (y_closest-y_pad)*(y_closest-y_pad) );
            //cout<<" PADPAIR: "<<x_closest<<" "<<y_closest<<" "<<x_pad<<" "<<y_pad<<" "<<dist<<endl;
            //cout<<j<<" "<<k<<endl;
            if(y_pad - y_closest > 0 ){
                distance = distance;         
            }else if(y_pad - y_closest < 0 ){
                distance = -distance; 
            }else{
                if(x_pad - x_closest > 0 ){
                    distance = distance;         
                }else if(x_pad - x_closest < 0 ){
                    distance = -distance; 
                } 
            }
            if (charge[j][k]!=0){
                //chi += ( PadPair[j][k].qfrac - func->Eval( PadPair[j][k].distance ) )*( PadPair[j][k].qfrac - func->Eval( PadPair[j][k].distance ) ) 
                //                                            / (   ( graph_FWs->Eval(PadPair[j][k].distance) )*( graph_FWs->Eval(PadPair[j][k].distance) )  );

                sigup = BMinvert( PadPair[j][k].qfrac + graph_FWs->Eval(distance)/2., distance );
                sigdn = BMinvert( PadPair[j][k].qfrac - graph_FWs->Eval(distance)/2., distance );
                exp_distance = BMinvert( PadPair[j][k].qfrac, distance );
                sig = graph_FWs->Eval(exp_distance)/func->Derivative(exp_distance);

                
                if( sig > 0.03 )  sig = 0.03;  //    sig = fabs(sigup-sigdn);   // sigma must no exceed 3 cm
                

                //cout<<distance<<"\t"<<exp_distance<<endl; 
                if (sig ==0) cout<<"############## SIGMA ERROR ##########"<<endl;


                chi += ( ( distance - exp_distance )  * ( distance - exp_distance ) ) 
                                                                                     / ( ( sig ) * ( sig ) );
                //cout<<"debug"<<endl;
                //cout<<x_pad<<" "<<y_pad<<endl;
                //cout<<" distance: "<<distance<<" qfrac meas: "<<PadPair[j][k].qfrac<<" exp dist: "<<exp_distance<<endl;
                //cout<<graph_FWs->Eval(distance)/2.<<"\t"<<sig;
            
                residual[j][k] = distance - exp_distance;
            }    
                //cout<<residual[j][k]<<endl;
            
        }
    }

    
            //cout<<"a,b,c: "<<a<<" "<<b<<" "<<c<<" chi: "<<chi<<endl;
            //cout<<PadPair[j][k].distance<<"\t"<<PadPair[j][k].qfrac<<"\t"<<sign<<endl;
            //cout << PadPair[j][k].distance << "\t" << BMinvert( PadPair[j][k].qfrac, PadPair[j][k].distance ) << endl;
            

    help->SetPoint(cnt,b,chi);
    std::ofstream of;
    of.open ("func.txt", std::ofstream::out | std::ofstream::app);
    of<<a<<" "<<b<<" "<<c<<" "<<chi<<endl;
    return chi;
    
}