/*
 * DCAP - Direct CAPture
 * Copyright (C) 2020 Oliver Kirsebom
 *
 * This program 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.
 *
 * This program 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 this program.  If not, see <https://www.gnu.org/licenses/>.
 */

#include <vector>
#include <ausa/constants/Mass.h>
#include <ORM/Particle.h>
#include <ORM/ParticleChannel.h>
#include <TF1.h>
#include <DCAP/DirectCapture.h>

using namespace std;
using namespace ORM;
using namespace DCAP;


int main(int argc, char* argv[]) {

    Particle p{1, 1, Parity::P, 0.5_J};   // proton
    Particle B11{5, 11, Parity::P, 1.5_J}; // B11 g.s. 3/2+

    double beamEnergy = 2000;
    
    vector<double> excEnergy{9640, 10844, 11820, 12710};
    vector<Spin> spin{3_J, 1_J, 2_J, 1_J};
    vector<OrbitalAngularMomentum> lOrb{2_L, 2_L, 2_L, 1_L};
    vector<double> S{0.3, 0.10, 0.15, 1}; //spectroscopic factors relative to 12.71
    for (auto& ss : S) ss *= 0.86;


    cout << Form("\nDC cross section @ %.2f MeV\n", beamEnergy/1e3) << endl;
    cout << " Ex (MeV)   l    S    cs (ub)" << endl;
    cout << "--------------------------------" << endl;

    for (size_t i=0; i<excEnergy.size(); i++) {

        vector<shared_ptr<ParticleChannel>> channels;
        for (size_t l=0; l<6; l++) {
            auto c = make_shared<ParticleChannel>(p, B11, spin[i], OrbitalAngularMomentum(l));
            c->setR0(1.36);
            channels.push_back(c);
        }

        auto& final = channels[lOrb[i].getValue()];

        double sigmaTot = 0;
        for (size_t l=0; l<6; l++) {
            double sigma = directCaptureCrossSec(channels[l], final, beamEnergy, excEnergy[i], S[i]);
            sigmaTot += sigma;
        } 

        cout << Form("  %.2f     %.0f   %.2f   %.2f", excEnergy[i]/1e3, lOrb[i].getValue(), S[i], sigmaTot) << endl;
    }
    cout << "--------------------------------\n" << endl;


    return 0;
}