/**
     @file NO_Thermal.cpp
     @author Lance Deaver

	$Date:$
    $Revision:$

    @copyright
     (©) Copyright 2009 by GATS Inc.
     11864 Canon Blvd., Suite 101, Newport News, VA 23606

     All Rights Reserved. No part of this software or publication may be
     reproduced, stored in a retrieval system, or transmitted, in any form
     or by any means, electronic, mechanical, photocopying, recording, or
     otherwise without the prior written permission of GATS Inc.
        
    @brief Corrects NO for the Thermal issues
     
*/

#include "NO_Thermal.h"
#include "NO_Correction.h"
#include "SigEarthLocParam.h"
#include "SOFIE_namespace.h"
#include "GATS_Utilities.hpp"

/*
#include <cmath>
#include <iostream>
#include <vector>
#include <stdexcept>
#include "SOFIE_namespace.h"
#include "GATS_Utilities.hpp"
*/

using namespace std;
using namespace GATS_Utilities;

NO_Thermal::NO_Thermal(Event& L0, ConfigFile& cf)
{


    if(cf.ValidEntry("SignalCorrection", "Max_ChiSquared") ) {
         maxChi =  cf.GetReal("SignalCorrection", "Max_ChiSquared");
         assert (maxChi >= 0. );
    } else {
       maxChi = 4.;
    }

	EventVar<double> startTime, endTime;
	L0.getEventVar("L0_BalanceStartTime",  startTime);
	L0.getEventVar("L0_BalanceEndTime", endTime);;

// since NO detector is Last we typically have a couple of extra seconds before NO starts
    BalanceStartTime = startTime[0] + 2.0;
    BalanceStopTime = endTime[0];


// std::cout << BalanceStartTime << "  " <<  BalanceStopTime << std::endl; exit(23);
}

//-----------------------------------------------------------------------------------------
//  Object method that sets up the NO thermal  correction for NO

//  Input:
//     inSignal........Input signals as multiple signals (All channels)
//
//  Output:
//     outSignal.......The corrected signals as multiple signals
//
//  Source: Lance Deaver, GATS
//-----------------------------------------------------------------------------------------
MultiSignal NO_Thermal::correctNO(MultiSignal &inSignal, SigEarthLocParam &earthInfo, double* resChi ) 
{
cout << "Inside NO_THermal" << endl;

//return inSignal;

    valarray<double> timeArray;
    vector<ChannelSignals> signalArray;
    Journal auditJournal= inSignal.getAudit();
    string comments;  
    string processName = "NO Thermal Correction";
    ChannelSignals Channel;


    try {
        inSignal.getTime(timeArray);    
	    vector<double> TP_Alts = earthInfo.getAltitudes();
	
	
		for(int i = 0; i< SOFIE::NO ; ++i) {
/*       for(int i = 0; i<= SOFIE::NO ; ++i) {

			Channel = inSignal.getChannel(i);
			valarray<double> strong = Channel.getSignal(SOFIE::Strong);   // get the strong band
            valarray<double> weak = Channel.getSignal(SOFIE::Weak);

            for(int j=0; j<strong.size(); ++j)
				if(timeArray[j] >= BalanceStartTime && timeArray[j] <= BalanceStopTime) 
					std::cout << timeArray[j] << "   " << strong[j] << std::endl;
std::cout << "&" << std::endl;
            for(int j=0; j<weak.size(); ++j)
                if(timeArray[j] >= BalanceStartTime && timeArray[j] <= BalanceStopTime)
                    std::cout << timeArray[j] << "   " << weak[j] << std::endl;
*/
		    signalArray.push_back( inSignal.getChannel(i) );
//std::cout << "&" << std::endl;
		}
//exit(23);	
	
        Channel = inSignal.getChannel(SOFIE::NO);   //NO channel
        valarray<double> strong = Channel.getSignal(SOFIE::Strong);   // get the strong band
		valarray<double> diffsig = Channel.getSignal(SOFIE::Diff); 

		int npts=(int)strong.size() ;  // number of points
        int iStart=0;  // index of starting location



// std::cout << *t_ptr << "  " << EventStopTime << "  " << timeArray[npts-1] << std::endl;
//std::cout << " npts " << npts << " size " << strong.size() << std::endl;
//exit(23);
//		int npts = strong.size();
		double maxBump= -0.0001 , weight=0.0 ; 
        double minFitz=140., maxFitz=190., startBA =  BalanceStartTime,
			endBA= BalanceStopTime; 
		char mode= (char)inSignal.getModeFlag(), band=(char)16;
	
		int chisFlag, bumpFlag;
	
		valarray<double> residual = strong , maxUncertianty(npts), 
	    	 minUncertianty(npts), Aparameters(9), Sigmaparameters(9) ;

// OK now call David's correction Code

//cout << " before davids code "   <<TP_Alts[0] <<  endl; exit(23);

//::sofieNO9_natural (
	::sofieNO12_natural (
		&strong[iStart], 
		&TP_Alts[iStart],
		&timeArray[iStart],
		&mode,
		&residual[iStart],
//		&band,
		&npts,
		&resChi[0],
		&maxChi,
		&chisFlag,
		&Aparameters[0],
		&Sigmaparameters[0],
		&maxUncertianty[0],
		&minUncertianty[0],
		&bumpFlag,
		&maxBump,
		&weight,
      &minFitz,
      &maxFitz,
      &startBA,
      &endBA );


//std::cout << resChi[0] << std::endl; exit(23);
if(chisFlag) {   // Bad fit
    cout << " **NO_Thermal Correction Not Performed**" << endl
     << "Chi Squared Value " << ConvertToString( *resChi) 
     << " Exceeds Limit: " << ConvertToString(maxChi) << endl;

 	cerr << " **NO_Thermal Correction Not Performed**" << endl 
     << "Chi Squared Value " << ConvertToString( *resChi) << " Exceeds Limit: " << ConvertToString(maxChi) ; 
	inSignal.getAudit().print(cerr);
	cerr << endl;

	resChi[0] *=  -1;

	return inSignal; 
}


/*
//std::cout << " before get sig " << std::endl;
valarray<double> r1 =  diffsig ; 
    ::sofieNO9_natural (
        &diffsig[iStart],
        &TP_Alts[iStart],
        &timeArray[iStart],
        &mode,
        &r1[iStart],
//      &band,
        &npts,
        &resChi[1],
        &maxChi,
        &chisFlag,
        &Aparameters[0],
        &Sigmaparameters[0],
        &maxUncertianty[0],
        &minUncertianty[0],
        &bumpFlag,
        &maxBump,
        &weight,
      &minFitz,
      &maxFitz,
      &startBA,
      &endBA );
if(chisFlag) {   // Bad fit
cout << " NO_Thermal correction Not Performed on Difference Signal " << endl;
	resChi[1] *= -1.;
	r1 = diffsig;
}
*/
//cout << "This is the number " <<  signalArray.size()  <<  " " << Channel.getID() << endl;

        signalArray.push_back( ChannelSignals(Channel.getID(), 
	             Channel.getAttenuation(SOFIE::Strong),
                 Channel.getAttenuation(SOFIE::Weak),
	             residual, 
                 Channel.getSignal(SOFIE::Weak),
	            //  r1 //  
                 Channel.getSignal(SOFIE::Diff)  
				)
	              );
		


	for(int i = SOFIE::NO+1 ; i< inSignal.numberOfSignals() ; ++i) {
	    signalArray.push_back( inSignal.getChannel(i) );
	}
		
/*
for(int i =0; i< residual.size(); ++i) 
    cout << ConvertToStringPrec(strong[i]) << "  " << ConvertToStringPrec(TP_Alts[i])  << endl;
cout << "&" << endl;

cout << "# " << ConvertToStringPrec( *resChi) << std::endl;
for(int i =0; i< residual.size(); ++i) 
    cout << ConvertToStringPrec(residual[i]) << "  " << ConvertToStringPrec(TP_Alts[i])  << endl;

exit(23);
 */   

    } catch (exception &ex) {
            
            cerr << "Exception thrown in NonLinearity :" << endl;
            cerr << "Channel Gas ID: " << Channel.getID() << endl;
            cerr << "Exception Message: " << ex.what() << endl;
            cerr << "Additional comments: " << comments << endl;
            cerr << "Current Audit Trail: ";
            inSignal.getAudit().print(cerr);
            cerr << endl;
            
            throw;
    }
    
     auditJournal.append(processName);
    
	cout << "finished NO_Thermal" << endl;
     return MultiSignal(inSignal.getEventNumber(), inSignal.get150kmTime(), 
                        inSignal.getModeFlag(), timeArray, signalArray, auditJournal) ;
}