//************************************************************************//
//									  //
//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)			  //
//  	      	   Julian Pychy (julian@ep1.rub.de)			  //
//          	   - Ruhr-Universit??t Bochum 				  //
//									  //
//  This file is part of Pawian.					  //
//									  //
//  Pawian is free software: you can redistribute it and/or modify	  //
//  it under the terms of the GNU General Public License as published by  //
//  the Free Software Foundation, either version 3 of the License, or 	  //
//  (at your option) any later version.	 	      	  	   	  //
//									  //
//  Pawian is distributed in the hope that it will be useful,		  //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of	  //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the	  //
//  GNU General Public License for more details.	      		  //
//									  //
//  You should have received a copy of the GNU General Public License     //
//  along with Pawian.  If not, see <http://www.gnu.org/licenses/>.	  //
//									  //
//************************************************************************//

// TensorPsiToGamXDecAmps class definition file. -*- C++ -*-
// Copyright 2013 Bertram Kopf

#include <getopt.h>
#include <fstream>
#include <string>

#include "PwaUtils/TensorPsiToGamXDecAmps.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "PwaUtils/DataUtils.hh"
#include "PwaUtils/IsobarTensorPsiToGamXDecay.hh"
#include "Particle/Particle.hh"
#include "ErrLogger/ErrLogger.hh"
#include "FitParams/AbsPawianParameters.hh"

TensorPsiToGamXDecAmps::TensorPsiToGamXDecAmps(std::shared_ptr<IsobarTensorPsiToGamXDecay> theDec, ChannelID channelID) :
  AbsXdecAmp(theDec, channelID)
  ,_noOfAmps(0)
{
  initialize();
  _noOfAmps=theDec->noOfAmplitudes();
  _ampLMap=theDec->ampLMap();

  _currentParamLocalMags.resize(_noOfAmps);
  _currentParamLocalPhis.resize(_noOfAmps);
  _currentParamLocalMagExpi.resize(_noOfAmps);
  _MagParamNames.resize(_noOfAmps);
  _PhiParamNames.resize(_noOfAmps);
  _MagParamNames[0]=_key+"Mag1";
  _PhiParamNames[0]=_key+"_1Phi";
  if(_noOfAmps>1){
    _MagParamNames[1]=_key+"Mag2";
    _PhiParamNames[1]=_key+"_2Phi";
  }
  if(_noOfAmps>2){
    _MagParamNames[2]=_key+"Mag3";
    _PhiParamNames[2]=_key+"_3Phi";
  }
}

TensorPsiToGamXDecAmps::TensorPsiToGamXDecAmps(std::shared_ptr<AbsDecay> theDec, ChannelID channelID) :
  AbsXdecAmp(theDec, channelID)
  ,_noOfAmps(0)
{
  initialize();
  //_noOfAmps=???
}

TensorPsiToGamXDecAmps::~TensorPsiToGamXDecAmps()
{
}


complex<double> TensorPsiToGamXDecAmps::XdecPartAmp(const Spin& lamX, Spin& lamDec, short fixDaughterNr, EvtData* theData, Spin& lamFs, AbsXdecAmp* grandmaAmp){
  
  complex<double> result(0.,0.);
  
  // Spin lam1Min=-_Jdaughter1;
  // Spin lam1Max= _Jdaughter1;
  // Spin lam2Min=-_Jdaughter2;
  // Spin lam2Max=_Jdaughter2;

  // if(fixDaughterNr == 1){
  //    lam1Min = lam1Max = lamDec;
  // }
  // else if(fixDaughterNr == 2){
  //    lam2Min = lam2Max = lamDec;
  // }
  // else{
  //    Alert << "Invalid fixDaughterNr in XdecPartAmp." << endmsg;
  // }

  // if(_enabledlamFsDaughter1){
  //   lam1Min=lamFs;
  //   lam1Max=lamFs;
  // }
  // else if(_enabledlamFsDaughter2){
  //   lam2Min=lamFs;
  //   lam2Max=lamFs;
  // }

  //  complex<double> result=lsLoop(grandmaAmp, lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);

  return result;
}




complex<double> TensorPsiToGamXDecAmps::XdecAmp(const Spin& lamX, EvtData* theData, AbsXdecAmp* grandmaAmp){
 
  complex<double> result(0.,0.);
  if( fabs(lamX) > _JPCPtr->J) return result;

  //  int evtNo=theData->evtNo;
  short currentSpinIndex=FunctionUtils::spin1IdIndex(_projId,lamX);
  
  if (!_recalculate){
    return _cachedAmpIdMap.at(theData->evtNo).at(_absDyn->grandMaId(grandmaAmp)).at(currentSpinIndex);
  }

  for(Spin lambda2=_lam2MinProj; lambda2<=_lam2MaxProj; ++lambda2){
    Id3StringType IdLamMotherLamGamLamX=FunctionUtils::spin3Index(lamX, _lam1MinProj, lambda2);
    complex<double> tmpResult(0.,0.);
    for(int i=0; i<_noOfAmps; ++i){
      tmpResult+=_currentParamLocalMagExpi.at(i)*theData->ComplexN3Spin.at(_decay->nameId()).at(i).at(IdLamMotherLamGamLamX)*_absDyn->eval(theData, grandmaAmp,_ampLMap.at(i));
    }

    result+=tmpResult*daughterAmp(lambda2, theData);    
  }

  if ( _cacheAmps){
     _cachedAmpIdMap[theData->evtNo][_absDyn->grandMaId(grandmaAmp)][currentSpinIndex]=result;
  }

  return result;
}

void  TensorPsiToGamXDecAmps::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){

  for (int i=0; i<_noOfAmps; ++i){
    //fill magnitude
    std::string magName=_MagParamNames.at(i);
    double valMag=1.;
    double errMag=0.5;
    //    double minMag=0.;
    //    double maxMag=30;

    fitPar->Add(magName, valMag, errMag);
    //    fitPar->SetLimits(magName, minMag, maxMag);

    //fill phi
    std::string phiName=_PhiParamNames.at(i);
    fitPar->Add(phiName, 0., 0.2);
  }

  _absDyn->fillDefaultParams(fitPar);
  if(!_daughter2IsStable) _decAmpDaughter2->fillDefaultParams(fitPar);
}

void TensorPsiToGamXDecAmps::fillParamNameList(){

  for (int i=0; i<_noOfAmps; ++i){
    //fill magnitude
    std::string magName=_MagParamNames.at(i);
    _paramNameList.push_back(magName);

    //fill phi
    std::string phiName=_PhiParamNames.at(i);
    _paramNameList.push_back(phiName);
  }
}

void TensorPsiToGamXDecAmps::print(std::ostream& os) const{
  return; //dummy
}


void TensorPsiToGamXDecAmps::updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar){

  for (int i=0; i<_noOfAmps; ++i){
    std::string magName=_MagParamNames.at(i);
    double theMag=std::abs(fitPar->Value(magName));
    
    std::string phiName=_PhiParamNames.at(i);
    double thePhi=fitPar->Value(phiName);
    _currentParamLocalMagExpi[i]=std::polar(theMag, thePhi);
  }

  _absDyn->updateFitParams(fitPar);
  if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(fitPar);
}

complex<double> TensorPsiToGamXDecAmps::daughterAmp(const Spin& lam2, EvtData* theData){
  complex<double> result(1.,0.);
  if(!_daughter2IsStable) result *= _decAmpDaughter2->XdecAmp(lam2, theData, this);
  return result;
}

