CascadeEvents.C

Go to the documentation of this file.

#ifndef __CINT__
#include <stdio.h>
#include <time.h>
#include "Riostream.h"
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <string>
#include <istream>
#include "TRandom.h"
#include "TTimeStamp.h"
#include "TMatrixD.h"
using namespace std;
#endif

using namespace Gw;

void cascataAGS(const Int_t, Int_t, FILE, const Double_t, const BaseGEM);

void matriceRES(FILE, Int_t*, Double_t**, Double_t*, Int_t, Int_t, Int_t*);

Int_t eccresiduoa(const Double_t&, Int_t&, Int_t&, Double_t&, const Double_t*);

Int_t eccresiduop(const Double_t&, Int_t&, Int_t&, Double_t&, const Double_t*);

Int_t eccresiduon(const Double_t&, Int_t&, Int_t&, Double_t&, const Double_t*);

Double_t trovame(Int_t&, Int_t&);

void posizioneAin_me();

void parabola(const Int_t*, const Double_t*, const Int_t&, Double_t&, Double_t&, Double_t&);

Double_t parabolame(Int_t&, Int_t&, Int_t&);

Double_t wei1(Int_t&, Int_t&);

void trovaAZ(const Int_t&, const Int_t&, Int_t&, Int_t&);


/*
//enrico
#define AV 15500.0                    
#define AS 16800.0
#define AC 720.0
#define APa 34000.0
#define AAS 23000.0
*/
#define MP 938280.0               
#define MN 939573.0
#define MALFA 3727409.0    

//frobrich
#define AV 15677.0                    
#define AS 18560.0
#define AC 720.0
#define APa 33500.0
#define AAS 23500.0


Int_t Z[3000], A[3000];
Double_t me[3000];
Int_t numnuclei;
Int_t posizA[3000];
Int_t AGS[200];
Int_t numA;
Int_t massa[29];
Int_t NmenoZ[11];
Int_t rettA, rettZ;
Int_t amin, zmin;


void CascadeEvents(const Char_t *finput, Int_t nbcas = 100, const Char_t *nameOUT = "toGEANT101.event")
{
        for (int i = 0; i < 3000; i++) {Z[i] = 0; A[i] = 0; me[i] = 0;}
        int i,j,k;
        for (i = 0; i < 200; i++) AGS[i] = 0;
        
        FILE *fileIN = fopen(finput,"r");
        if ( fileIN == NULL ) {
                printf(" Cannot open file %s \n",finput); 
                return;
        }

        //---------carico le masse
        FILE *fileMAS = fopen("mastab.dat","r");
        if ( fileMAS == NULL ) {
                printf(" Cannot open file mastab.dat \n"); 
                return;
        }
        Char_t Line[300];
        fgets (Line, 300, fileMAS);
        sscanf(Line,"%*c%*c%*d%*d%d%d%*s%*s%lf", &Z[0], &A[0], &me[0]);
        i = 0;
        while (!feof(fileMAS)){
                fgets (Line, 300, fileMAS);
                sscanf(Line,"%*c%*c%*d%*d%d%d%*s%lf", &Z[i], &A[i], &me[i]);
                if (A[i] > 0) i++;
        }
        numnuclei = i;
        posizioneAin_me();

        //---------carico gli spettri evaporativi
        Char_t buf[100], buf2[100], buf3[100];
        while ((strcmp(buf, "SUMMED")) || (strcmp(buf2, "EVAPORATION")) || (strcmp(buf3, "SPECTRA"))) {
                fgets (Line, 300, fileIN);
                sscanf(Line,"%s %s %s", buf, buf2, buf3);
        }
        for (int i=0; i < 3; i++) fgets (Line, 300, fileIN);
        Float_t epart[50], dn[50], dp[50], da[50], eg[50], dg[50];
        int indice=0;
        while (strcmp(buf, "0")) {
                fgets (Line, 300, fileIN);
                sscanf(Line,"%s", buf);
                if (strcmp(buf, "0")) {
                        sscanf(Line,"%f%f%f%f %*41c %f %f", &epart[indice], &dn[indice], &dp[indice], &da[indice], &eg[indice], &dg[indice]);
                }
                indice++;
        }
        Float_t dcumn[50], dcump[50], dcuma[50], dcumg[50];
        dcumn[0] = 0; dcump[0] = 0; dcuma[0] = 0; dcumg[0] = 0;
        for (int i = 1; i < 50; i++){
                dcumn[i] = dcumn[i-1] + dn[i-1];
                dcump[i] = dcump[i-1] + dp[i-1];
                dcuma[i] = dcuma[i-1] + da[i-1];
        }
        
        //---------calcolo l'enegia "media" delle particelle evaporate (en media := en t.c. il 50% ha meno energia)
        Float_t eamed, epmed, enmed;
        Float_t soglia = dcuma[49]*0.5;
        for (int j = 0; j < 50; j++){
                if (dcuma[j] < soglia) eamed = 1000 * epart[j];
        }
        soglia = dcump[49]*0.5;
        for (int j = 0; j < 50; j++){
                if (dcump[j] < soglia) epmed = 1000 * epart[j];
        }
        soglia = dcumn[49]*0.5;
        for (int j = 0; j < 50; j++){
                if (dcumn[j] < soglia) enmed = 1000 * epart[j];
        }
        cout << "eamed = " << eamed << " epmed = " << epmed << " enmed = " << enmed << '\n';
        
        //---------carico la reazione
        while (strcmp(buf, "ENTR.CHANNEL,")) {
                fgets (Line, 300, fileIN);
                sscanf(Line,"%s", buf);
        }
        Int_t zp, ap, zb, ab, zcn, acn;
        Float_t ep, ecccn;
        sscanf(Line,"%*s%*s%*s %d%d %d%d %f %d%d %f", &zp, &ap, &zb, &ab, &ep, &zcn, &acn, &ecccn);

        //---------carico la matrice dei nuclei residui
        //Int_t massa[29]; //globale
        Double_t pop[29][6];
        Double_t popsum[29];
        Int_t numpari, numdispari;
        //Int_t NmenoZ[11]; //globale
        matriceRES(fileIN, massa, pop, popsum, numpari, numdispari, NmenoZ);  //carica la matrice con i nuclei residui

        Double_t total = 0;
        Int_t numpop=0;
        for (i = 0; i < 29; i++){
                if (massa[i] % 2) {for (j = 0; j < numdispari; j++) {total += pop[i][j]; if (pop[i][j] > 0) numpop++; }} //se massa[i] e' dispari
                else {for (j = 0; j < numpari; j++) {total += pop[i][j]; if (pop[i][j] > 0) numpop++;}} //se massa[i] e' pari
        }

        Double_t popcum[29][6];
        Double_t sumparz = 0;
        for (i = 0; i < 29; i++) {
                if (massa[i] % 2) {for (j = 0; j < numdispari; j++) {sumparz += pop[i][j] / total; popcum[i][j] = sumparz;}} //se massa[i] e' dispari
                else {for (j = 0; j < numpari; j++) {sumparz += pop[i][j] / total; popcum[i][j] = sumparz;}} //se massa[i] e' pari
        }
        
        //---------trovo amin e zmin dei nuclei popolati
        Int_t atmp, ztmp;
        amin = 99999;
        zmin = 99999;
        for (i = 0; i < 29; i++){
                if (massa[i] % 2) {for (j = 0; j < numdispari; j++) {if (pop[i][j] > 0) { trovaAZ(i, j, atmp, ztmp); if (atmp < amin) amin = atmp; if (ztmp < zmin) zmin = ztmp;}}}//se massa[i] e' dispari
                else {for (j = 0; j < numpari; j++) {if (pop[i][j] > 0) { trovaAZ(i, j, atmp, ztmp); if (atmp < amin) amin = atmp; if (ztmp < zmin) zmin = ztmp;}}} //se massa[i] e' pari
        }

        //---------trovo gli eccessi di massa nel rettangolo (amin,zmin)-(acn,zcn)
        rettA = acn-amin+1;
        rettZ = zcn-zmin+1;
        Double_t *mecalc = new Double_t[rettA * rettZ];
        cout << "Calcolo la massa dei nuclei coinvolti ...\n";
        for (int a = 0; a < rettA; a++) {
                for (int z = 0; z < rettZ; z++) {
                        mecalc[a * rettZ + z] = trovame(a + amin, z + zmin);
                        //cout << "mecalc[" << (a + amin) << "][" << (z + zmin) << "]: " << mecalc[a * rettZ + z] << '\n';
                }
        }
        cout << "Fatto\n";
        
        //---------inizializzo gli spettri
        cout << "spettri\n";
        Int_t r=0;
        Int_t residuo[165];
        for (i = 0; i < 29; i++){
                if (massa[i] % 2) {for (j = 0; j < numdispari; j++) {if (pop[i][j] > 0) {residuo[r] = i * 10 + j; r++;}}} //se massa[i] e' dispari
                else {for (j = 0; j < numpari; j++) {if (pop[i][j] > 0) {residuo[r] = i * 10 + j; r++;}}} //se massa[i] e' pari
        }
        for (int enr=0; enr<r; enr++) {cout << enr << " - " << residuo[enr] << '\n';}
        cout << "spettri\n";
        
        // load AGS file(s)
        BaseGEM* random = new BaseGEM[numpop];
        cout << "numpop: " << numpop << "\n";
        for (int r = 0; r < numpop; r++) {
                //trovo A e Z del nucleo r-esimo
                Int_t posi = residuo[r] / 10;
                Int_t posj = residuo[r] - 10 * posi;
                Int_t posA, posZ;
                trovaAZ(posi, posj, posA, posZ);
                /*
                Int_t posA = massa[posi];
                Int_t posNmZ;
                if ((massa[posi] % 2) == (NmenoZ[0] % 2)) {posNmZ = NmenoZ[2 * posj];}
                else {posNmZ = NmenoZ[2 * posj + 1];}
                Int_t posZ = (posA - posNmZ) / 2;
                */
                
                //converto nell'elemento
                //char chem[109][3] = {"H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt"};
                char chem[259][3] = {"H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X","X"};

                //costruisco il nomefileAGS
                char nomefileAGS [500];
                sprintf (nomefileAGS, "../ags/%d%s_ensdf_A.ags", posA, chem[posZ-1]);
                printf ("nomefileAGS: %s\n", nomefileAGS);
                Char_t *nameAGS = nomefileAGS;

                //controllo se il file AGS esiste               
                FILE *fileAGS = fopen(nameAGS,"r");
                if ( fileAGS == NULL ) {
                        printf(" File AGS non esiste\n"); 
                        AGS[r] = 0;
                }
                else {                          
                        if ( random[r].InitAGS(nameAGS) == 0 ) {
                                cout << " GLS level Scheme not correctly loaded " << endl << "Nucleo Z:" << posZ << " A:" << posA << " posi:" << posi << " posj:" << posj << endl; 
                                AGS[r] = 0;
                        }
                        else {
                                cout << "Loaded GLS level Scheme for nucleo " << r << '\n';
                                AGS[r] = 1;
                        }
                }
                cout << posA << chem[posZ-1] << " (" << posZ << ") in posizione " << residuo[r] << " processato" << '\n';
                cout << "numpop: " << numpop << " |r: " << r << '\n';
                if (r == (numpop -1)) break;
        }
        
        cout << numpop << " FILE AGS CARICATI" << '\n';

        
        // open a file to store the simulated events.
        FILE *fileOUT = fopen(nameOUT,"w");
        if ( fileOUT == NULL ) {
                printf(" Cannot open file %s \n",nameOUT); return;
        }
        fprintf(fileOUT,"FORMAT 4 4 \n");
        fprintf(fileOUT,"# \n");
        fprintf(fileOUT,"# PARTICLES GENERATED FROM CASCADE OUTPUT \n");
        fprintf(fileOUT,"# GAMMAS GENERATED FROM THE LEVEL SCHEME \n");
        fprintf(fileOUT,"# \n");
        fprintf(fileOUT,"REACTION %d %d %d %d %.1f\n", zp, ap, zb, ab, ep);
        fprintf(fileOUT,"# \n");
        fprintf(fileOUT,"EMITTED 4 1 2 3 7 \n");
        fprintf(fileOUT,"# \n");
        
        
        
        
        
        //---------genero eventi
        Int_t numev[29][6];
        //inizializzo gRandom con i nanosecondi
        Double_t rndmdb, rndmener;
        Double_t nanos = TTimeStamp().GetNanoSec();
        //gRandom -> SetSeed((UInt_t)nanos);
        gRandom -> SetSeed(0);
        for (i = 0; i < 29; i++) for (j = 0; j < 6; j++) numev[i][j] = 0;
        //Float_t eaemitted[10], epemitted[10], enemitted[10];
        Int_t Aemitter,Zemitter;
        Double_t eccemitter;

        //******************************************************************************************************
        //*******************************INIZIO***********************************************************************
        //******************************************************************************************************

        for (k = 0; k < nbcas; k++) {
                //------------------------estraggo un nucleo residuo
                rndmdb = gRandom -> Rndm();
                Int_t posi = 0;
                Int_t posj = 0;
                Int_t posOK = 0;
                for (i = 0; i < 29; i++) {
                        if (massa[i] % 2) {for (j = 0; j < numdispari; j++) {if ((rndmdb < popcum[i][j]) && (posOK == 0)) {posOK = 1; posi = i; posj = j;}}} //se massa[i] e' dispari
                        else {for (j = 0; j < numpari; j++) {if ((rndmdb < popcum[i][j]) && (posOK == 0)) {posOK = 1; posi = i; posj = j;}}} //se massa[i] e' pari
                        if (posOK) break;
                }


                numev[posi][posj]++;
                Int_t posA, posZ;
                trovaAZ(posi, posj, posA, posZ);
                /*
                //trovo A
                Int_t posA;
                posA = massa[posi];
                //trovo N-Z
                Int_t posNmZ;
                if ((massa[posi] % 2) == (NmenoZ[0] % 2)) {posNmZ = NmenoZ[2 * posj];}
                else {posNmZ = NmenoZ[2 * posj + 1];}
                Int_t posZ = (posA - posNmZ) / 2;
                */
                
                //------------------------numero di particelle da emettere
                Int_t diffp = zcn - posZ;
                Int_t diffn = acn - posA - diffp;
                Int_t nalfa = 0, nprot = 0, nneut = 0;
                while ((diffp >= 2) && (diffn >= 2)) {
                        nalfa++;
                        diffp -= 2;
                        diffn -= 2;
                }
                nprot = diffp;
                nneut = diffn;
                //cout << "Acn: " << acn << " Zcn: " << zcn << " Af: " << posA << " Zf: " << posZ << " alfa: " << nalfa << " prot: " << nprot << " neut: " << nneut << '\n';
                
                
                //------------------------emetti particelle
                Double_t* eaemitted = new Double_t[nalfa];  //Double_t *eaemitted; //eaemitted = (Double_t*)calloc(nalfa, sizeof(Double_t));
                Double_t* epemitted = new Double_t[nprot];
                Double_t* enemitted = new Double_t[nneut];
                Int_t ok = 1;
                Int_t risparopart = 0;
                Double_t ridprob = 1;
                
                while (ok == 1){
                        ok = 0;
                        Aemitter = acn;
                        Zemitter = zcn;
                        eccemitter = ecccn * 1000;
                        //cout << "eccemitter = " << eccemitter << '\n';
                        for (int i = 0; i < nalfa; i++){
                                if (eccemitter <= eamed) ok = 1;
                                //rndmener = dcuma[49]*(gRandom -> Rndm());
                                rndmener = dcuma[49]*(gRandom -> Rndm()) * ridprob;
                                for (int j = 0; j < 50; j++){
                                        if (dcuma[j] > rndmener) {
                                                eaemitted[i] = 1000 * (epart[j-1] - 0.5 + (rndmener - dcuma[j-1]) / da[j-1]);
                                                //cout << "mecalc " << mecalc << "  &mecalc " << &mecalc << "\n";
                                                eccresiduoa(eaemitted[i], Aemitter, Zemitter, eccemitter, mecalc);
                                                break;
                                        }
                                }
                                //cout << "eccemitter = " << eccemitter << " eaemitted = " << eaemitted[i] << '\n';
                        }
                        for (int i = 0; i < nprot; i++){
                                if (eccemitter <= epmed) ok = 1;
                                //rndmener = dcump[49]*(gRandom -> Rndm());
                                rndmener = dcump[49]*(gRandom -> Rndm()) * ridprob;
                                for (int j = 0; j < 50; j++){
                                        if (dcump[j] > rndmener) {
                                                epemitted[i] = 1000 * (epart[j-1]  - 0.5 + (rndmener - dcump[j-1]) / dp[j-1]);
                                                //cout << "mecalc " << mecalc << "  &mecalc " << &mecalc << "\n";
                                                eccresiduop(epemitted[i],Aemitter,Zemitter,eccemitter, mecalc);
                                                break;
                                        }
                                }
                                //cout << "eccemitter = " << eccemitter << " epemitted = " << epemitted[i] << '\n';
                        }
                        for (int i = 0; i < nneut; i++){
                                if (eccemitter <= enmed) ok = 1;
                                //rndmener = dcumn[49]*(gRandom -> Rndm());
                                rndmener = dcumn[49]*(gRandom -> Rndm()) * ridprob;
                                for (int j = 0; j < 50; j++){
                                        if (dcumn[j] > rndmener) {
                                                enemitted[i] = 1000 * (epart[j-1]  - 0.5 + (rndmener - dcumn[j-1]) / dn[j-1]);
                                                //cout << "mecalc " << mecalc << "  &mecalc " << &mecalc << "\n";
                                                eccresiduon(enemitted[i],Aemitter,Zemitter,eccemitter, mecalc);
                                                break;
                                        }
                                }
                                //cout << "eccemitter = " << eccemitter << " enemitted = " << enemitted[i] << '\n';
                        }
                        if (eccemitter <= 0) ok = 1;
                        if (ok == 1) {risparopart++; ridprob *= 0.9; /*cout <<"RISPARO PARTICELLE (energia)************" << risparopart << '\n';*/}
                }
                
                //------------------------scrivi dati sorgente
                fprintf(fileOUT," $ \n -101 \n ");
                
                //------------------------scrivi le particelle
                for (int i = 0; i < nalfa; i++) fprintf(fileOUT,"7 %.1lf 0. 0.\n ", eaemitted[i]);
                for (int i = 0; i < nprot; i++) fprintf(fileOUT,"3 %.1lf 0. 0.\n ", epemitted[i]);
                for (int i = 0; i < nneut; i++) fprintf(fileOUT,"2 %.1lf 0. 0.\n ", enemitted[i]);
                
                /*//converto nell'elemento
                char chem[109][3] = {"H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt"};
                //costruisco il nomefileAGS
                char nomefileAGS [500];
                sprintf (nomefileAGS, "../ags/%d%s_ensdf_A.ags", posA, chem[posZ-1]);
                printf ("nomefileAGS: %s", nomefileAGS);
                //sparo una cascata
                Char_t *nameAGS = nomefileAGS;*/
                
                //------------------------emetti e scrivi gamma nello spettro noto
                //trovo r(posi,posj) -> resi
                Int_t resi = 0;
                while (residuo[resi] != (10 *posi + posj)) resi++;
                cascataAGS(resi, 1, fileOUT, eccemitter, random);

                if(k && !(k%1000)) {printf("%10d eventi\r",k);}
        }

        //******************************************************************************************************
        //*******************************FINE***********************************************************************
        //******************************************************************************************************

        
        //output matrice (N-Z) - A dei nuclei simulati (frequenze)
        Int_t jmax;
        for (i = 0; i < 29; i++) {
                if ((massa[i] % 2) == ((NmenoZ[0] * NmenoZ[0]) % 2)) {jmax = 6;} 
                else {jmax = 5;} 
                for (j = 0; j < jmax; j++) {cout << numev[i][j] << "\t";}
                cout << '\n';
        }

        // close the output file
        fprintf(fileOUT, "# \n"); 
        fclose(fileOUT);
        cout << '\n';
}




void cascataAGS(const Int_t r, Int_t nbcas = 1, FILE *fileOUT, const Double_t eccemitter, const BaseGEM random[])
{

        //cout << "nucleo " << r << '\n';
        // start the Monte-Carlo
        Cascade cas; GammaLink *gam;
        Double_t etotg;
        Int_t risparogamma = 0;
        for (Int_t i = 0; i < nbcas; i++ ) {
                if ((AGS[r] == 1) && (risparogamma < 50)) {
                        etotg = 0;
                        random[r].DoCascade(&cas); // next cascade
                        for (Int_t j = 0; j < cas.GetSize(); j++ ) {
                                gam = (GammaLink *)cas.At(j); 
                                etotg += gam->GetEnergy().Get();
                        }
                        //cout << "etotg = " << etotg << " |eccemitter = " << eccemitter << '\n';
                        if (etotg >= eccemitter) {
                                //cout << "RISPARO GAMMA" << '\n';
                                risparogamma++;
                                i--;
                                continue;
                        }
                        else {
                                risparogamma = 0;
                                fprintf(fileOUT,"1 %.1f 0. %.1f \n ", (eccemitter - etotg), 0); //gamma di raccordo
                                for (Int_t j = 0; j < cas.GetSize(); j++ ) {
                                        gam = (GammaLink *)cas.At(j); 
                                        fprintf(fileOUT,"1 %.1lf 0. %.1f \n ",gam->GetEnergy().Get(),0.);
                                }
                                //cout << "gammascritti" << "\n\n";
                        }
                }
                else {
                        fprintf(fileOUT,"1 %.1lf 0. %.1f \n ", eccemitter, 0.);//gamma di raccordo e basta
                        //cout << "gammascritto" << "\n\n";
                }
        }
}



void matriceRES(FILE *fileIN, Int_t massa[29], Double_t pop[29][6], Double_t popsum[29], Int_t& numpari, Int_t& numdispari, Int_t NmenoZ[11])
{
        //mi porto su OUTPUT OPTIONS
        char s[256];
        //int num;
        fscanf(fileIN, "%s", s); 
        Int_t outputopt = 0;
        while (outputopt == 0) {
                while (strcmp(s, "OUTPUT")) fscanf(fileIN, "%s", s);
                fscanf(fileIN, "%s", s);
                if (!strcmp(s, "OPTIONS")) {outputopt = 1;}
        }
        
        //mi porto su ------...
        for (int i=0; i < 15; i++) {
                fscanf(fileIN, "%s",s);
                if (!strcmp(s, "------------------------------------------------------------------------------------")) {break;}
        }

        fscanf(fileIN, "%s", s); // carico A
        //Int_t NmenoZ[11];
        for (int i = 0; i < 11; i++) {fscanf(fileIN, "%d", &NmenoZ[i]);}
        //Int_t numpari, numdispari;
        
        if ((NmenoZ[0] * NmenoZ[0]) % 2) { // se NmenoZ[0] e' dispari
                numpari = 5;    
                numdispari = 6; 
        }
        else {
                numpari = 6;    
                numdispari = 5; 
        }
        /*
        fscanf(fileIN, "%s", &s); // carico "(N-Z)"
                
        for (i = 0; i < 29; i++) {
                fscanf(fileIN, "%d", &massa[i]); 
                if (massa[i] % 2) {for (int j = 0; j < numdispari; j++) {fscanf(fileIN, "%lf", &pop[i][j]);}} //se massa[i] e' dispari
                else {for (j = 0; j < numpari; j++) {fscanf(fileIN, "%lf", &pop[i][j]);}} //se massa[i] e' pari
                fscanf(fileIN, "%lf", &popsum[i]); 
                fscanf(fileIN, "%s", &s); //carico il grafico
        }
        */

        Char_t Line[300];
        fgets (Line, 300, fileIN);
        fgets (Line, 300, fileIN);
        for (int i = 0; i < 29; i++) {
                fgets (Line, 300, fileIN);
                sscanf(Line,"%d",&massa[i]);
                if ((massa[i]%2) == ((NmenoZ[0]*NmenoZ[0])%2)) sscanf(Line,"%*d%lf%lf%lf%lf%lf%lf%lf", &pop[i][0], &pop[i][1], &pop[i][2], &pop[i][3], &pop[i][4], &pop[i][5], &popsum[i]);
                else sscanf(Line,"%*d%lf%lf%lf%lf%lf%lf", &pop[i][0], &pop[i][1], &pop[i][2], &pop[i][3], &pop[i][4], &popsum[i]);
        }
        //cout << "esco da matriceres\n";               
        
}



Int_t eccresiduoa(const Double_t &eaemitted, Int_t &Aemitter, Int_t &Zemitter, Double_t &eccemitter, const Double_t* mecalc)
{
        //cout << "mecalc " << mecalc << "  *mecalc " << *mecalc << "\n\n";
        Double_t u = 931494.043;
        Double_t ma = 3728401.08760000021;
        //Double_t meemitter = trovame(Aemitter, Zemitter);
        Double_t meemitter = mecalc[(Aemitter - amin) * rettZ + (Zemitter - zmin)];
        //cout << "A: " << Aemitter << "  Z: " << Zemitter << "massa emitter: " << (meemitter + Aemitter * u)<< '\n';
        //Double_t meresiduo = trovame((Aemitter - 4), (Zemitter - 2));
        Double_t meresiduo = mecalc[(Aemitter - 4 - amin) * rettZ + (Zemitter - 2 - zmin)];
        //cout << "A: " << (Aemitter - 4) << "  Z: " << (Zemitter - 2) << "massa residuo: " << (meresiduo + (Aemitter - 4) * u) << '\n';
        eccemitter = eccemitter + Aemitter * u + meemitter - eaemitted - ma - sqrt(((Aemitter-4)*u + meresiduo) * ((Aemitter-4)*u + meresiduo) + eaemitted * eaemitted + 2 * eaemitted * ma);
        Aemitter -= 4;
        Zemitter -= 2;
        return 0;
}



Int_t eccresiduop(const Double_t &epemitted, Int_t &Aemitter, Int_t &Zemitter, Double_t &eccemitter, const Double_t* mecalc)
{
        //cout << "mecalc " << mecalc << "  *mecalc " << *mecalc << "\n\n";
        Double_t u = 931494.043;
        Double_t mp = 938783.013500000038;
        //Double_t meemitter = trovame(Aemitter, Zemitter);
        Double_t meemitter = mecalc[(Aemitter - amin) * rettZ + (Zemitter - zmin)];
        //cout << "A: " << Aemitter << "  Z: " << Zemitter << "massa emitter: " << (meemitter + Aemitter * u)<< '\n';
        //Double_t meresiduo = trovame((Aemitter - 1), (Zemitter - 1));
        Double_t meresiduo = mecalc[(Aemitter - 1 - amin) * rettZ + (Zemitter - 1 - zmin)];
        //cout << "A: " << (Aemitter - 1) << "  Z: " << (Zemitter - 1)<< "massa residuo: " << (meresiduo + (Aemitter - 1) * u) << '\n';
        eccemitter = eccemitter + Aemitter * u + meemitter - epemitted - mp - sqrt(((Aemitter-1)*u + meresiduo) * ((Aemitter-1)*u + meresiduo) + epemitted * epemitted + 2 * epemitted * mp);
        Aemitter -= 1;
        Zemitter -= 1;
        return 0;
}



Int_t eccresiduon(const Double_t &enemitted, Int_t &Aemitter, Int_t &Zemitter, Double_t &eccemitter, const Double_t* mecalc)
{
        //cout << "mecalc " << mecalc << "  *mecalc " << *mecalc << "\n\n";
        Double_t u = 931494.043;
        Double_t mn = 939565.360100000021;
        //Double_t meemitter = trovame(Aemitter, Zemitter); //keV
        Double_t meemitter = mecalc[(Aemitter - amin) * rettZ + (Zemitter - zmin)];
        //cout << "A: " << Aemitter << "  Z: " << Zemitter << "massa emitter: " << (meemitter + Aemitter * u)<< '\n';
        //Double_t meresiduo = trovame((Aemitter - 1), Zemitter);
        Double_t meresiduo = mecalc[(Aemitter - 1 - amin) * rettZ + (Zemitter - zmin)];
        //cout << "A: " << (Aemitter - 1) << "  Z: " << Zemitter << "massa residuo: " << (meresiduo + (Aemitter - 1) * u) << '\n';
        eccemitter = eccemitter + Aemitter * u + meemitter - enemitted - mn - sqrt(((Aemitter-1)*u + meresiduo) * ((Aemitter-1)*u + meresiduo) + enemitted * enemitted + 2 * enemitted * mn);
        Aemitter -= 1;
        return 0;
}



Double_t trovame(Int_t &Anuc, Int_t &Znuc)
{ 
        Double_t massapar;
        //me tabulato
        for (int z = posizA[Anuc]; z < numnuclei; z++) {
                if ((Znuc == Z[z]) && (Anuc == A[z])) {/*cout << "me tabulato     ";*/ return me[z];}
        }
        
        //me non tabulato ma ho nuclei con lo stesso A
        Int_t numA = 0;
        for (int z = posizA[Anuc]; z < numnuclei; z++) {
                if (A[z] == Anuc) numA++; //conta nuclei con massa Anuc
                else break;
        }
        //cout << "posizA Anuc " << posizA[Anuc]<< "    numA " << numA << '\n';
        if (numA >= 3) {
                massapar = parabolame(Anuc, Znuc, numA);
                //cout << "massapar " << Znuc << " " << massapar << '\n';
                //cout << "me da parabola  ";
                return massapar;
        }
        
        //me da stimare
        Double_t massawei;
        massawei = wei1(Anuc, Znuc);
        //cout << "me da weizsacker  ";
        return massawei;
        
}



void posizioneAin_me()
{
        cout << "numnuclei " << numnuclei << " | ";
        for (int a = 1; a < A[numnuclei - 1]; a++) {
                for (int i = 0; i < numnuclei; i++) {
                        if (A[i] == a) {posizA[a] = i; break;}
                }
        }
        cout << "posizioneAin_me OK" << '\n';
}



void parabola(const Int_t x[], const Double_t y[], const Int_t &Nnuc, Double_t &a, Double_t &b, Double_t &c)
{
        Int_t* x2 = new Int_t[Nnuc];
        Int_t* x3 = new Int_t[Nnuc];
        Int_t* x4 = new Int_t[Nnuc];
        Double_t* xy = new Double_t[Nnuc];
        Double_t* x2y = new Double_t[Nnuc];
        for (int i = 0; i < Nnuc; i++) {
                x2[i] = x[i] * x[i];
                x3[i] = x2[i] * x[i];
                x4[i] = x3[i] * x[i];
                xy[i] = x[i] * y[i];
                x2y[i] = xy[i] * x[i];
        }
        Int_t Sx, Sx2, Sx3, Sx4;
        Double_t Sy, Sxy, Sx2y;
        for (int i = 0; i < Nnuc; i++) {
                Sx += x[i];
                Sx2 += x2[i];
                Sx3 += x3[i];
                Sx4 += x4[i];
                Sy += y[i];
                Sxy += xy[i];
                Sx2y += x2y[i];
        }
        Double_t sumel[9];
        sumel[0] = Sx4;
        sumel[1] = Sx3;
        sumel[2] = Sx2;
        sumel[3] = Sx3;
        sumel[4] = Sx2;
        sumel[5] = Sx;
        sumel[6] = Sx2;
        sumel[7] = Sx;
        sumel[8] = Nnuc;
        TMatrixD sum(3, 3, sumel,0);
        /*cout << sum[0][0] << " " << sum[0][1] << " " << sum[0][2] << "\n";
        cout << sum[1][0] << " " << sum[1][1] << " " << sum[1][2] << "\n";
        cout << sum[2][0] << " " << sum[2][1] << " " << sum[2][2] << "\n";*/
        Double_t det;
        TMatrixD sumold = sum;
        sum.Invert(&det);
        
        /*
        cout << "\n MATRICI\n";
        cout << sum[0][0] << " \t" << sum[0][1] << " \t" << sum[0][2] <<'\n';
        cout << sum[1][0] << " \t" << sum[1][1] << " \t" << sum[1][2] <<'\n';
        cout << sum[2][0] << " \t" << sum[2][1] << " \t" << sum[2][2] <<"\n\n";
        
        cout << sumold[0][0] << " \t" << sumold[0][1] << " \t" << sumold[0][2] <<'\n';
        cout << sumold[1][0] << " \t" << sumold[1][1] << " \t" << sumold[1][2] <<'\n';
        cout << sumold[2][0] << " \t" << sumold[2][1] << " \t" << sumold[2][2] <<"\n\n";
        
        TMatrixD prod(3, 3);
        for (int h = 0; h < 3; h++){
                for (int k = 0; k < 3; k++){
                        for (int i = 0; i < 3; i++){
                                prod[h][k] += sum[h][i] * sumold[i][k];  
                        }
                }
        }
        cout << prod[0][0] << " \t" << prod[0][1] << " \t" << prod[0][2] <<'\n';
        cout << prod[1][0] << " \t" << prod[1][1] << " \t" << prod[1][2] <<'\n';
        cout << prod[2][0] << " \t" << prod[2][1] << " \t" << prod[2][2] <<'\n';
        */
        
        a = sum[0][0] * Sx2y + sum[0][1] * Sxy + sum[0][2] * Sy;
        b = sum[1][0] * Sx2y + sum[1][1] * Sxy + sum[1][2] * Sy;
        c = sum[2][0] * Sx2y + sum[2][1] * Sxy + sum[2][2] * Sy;
}


Double_t parabolame(Int_t &Anuc, Int_t &Znuc, Int_t &numA)
{
        Double_t a, b, c, massa;

        switch (Anuc%2) {
        case 0: //pari-pari e diapari-dispari
                int numAp = 0;
                Int_t* Znp = new Int_t[numA];
                Double_t* menp = new Double_t[numA];
                
                for (int i = 0; i < numA; i++) {        
                        if ((Z[posizA[Anuc] + i]%2) == Znuc%2) {
                                Znp[numAp] = Z[posizA[Anuc] + i]; 
                                menp[numAp] = me[posizA[Anuc] + i];
                                numAp++;
                        }
                }
                parabola(Znp, menp, numAp, a, b, c);
                massa = a * Znuc * Znuc + b * Znuc + c;
                //cout << "PARI-PARI DISP-DISP" << "PARI-PARI DISP-DISP" << "PARI-PARI DISP-DISP" << '\n'; 
                return massa;
                
        case 1: //pari-dispari
                Int_t* Zn = new Int_t[numA];
                Double_t* men = new Double_t[numA];
                //cout << "numA " << numA<< '\n';
                //for (int i = 0; i < numA; i++) cout << "Z[posizA[Anuc] + i]" << Z[posizA[Anuc] + i] << '\n';
                for (int i=0; i<numA; i++){
                        Zn[i] = Z[posizA[Anuc] + i];
                        men[i] = me[posizA[Anuc] + i];
                }
                //for (int i=0; i<numA; i++) cout << Zn[i] << " " << men[i] << "|";
                parabola(Zn, men, numA, a, b, c);
                massa = a * Znuc * Znuc + b * Znuc +c;
                //cout << "PARI-DISP"  << "PARI-DISP" << "PARI-DISP" << "PARI-DISP" << "PARI-DISP" << "PARI-DISP" << '\n';
                return massa;
        }
        
}


Double_t wei1(Int_t &a, Int_t &z)
{
  double mas;
  double delt;
  
  mas = AV*a;
  mas -= AS * pow(a, (2./3.));
  mas -= AC * z * (z-1)/pow(a, (1./3.));
  mas -= AAS * pow((a-2*z), 2.)/a;
  /*
  if( (a/2-floor(a/2)) > 0.30 ) delt = 0.;
  else if( (z/2-floor(z/2)) > 0.30 ) delt = -APa * pow(a, -0.75);
        else delt =  APa * pow(a, -0.75);
  */
        if (a%2 == 1) delt = 0.;
        else if (z%2 == 1) delt = -APa * pow(a, -0.75);
                else delt =  APa * pow(a, -0.75);
                
  mas += delt;
  mas -= 2*mas;
  mas += z*MP;
  mas += (a-z)*MN;
  
        mas -= a * 931494.043;
  return mas;
}


void trovaAZ(const Int_t &posi, const Int_t &posj, Int_t &posA, Int_t &posZ)
{
        //for (int i = 0; i < 11; i++) cout << "NmenoZ[" << i << "] = " << NmenoZ[i] << '\n'; 
        //cout << "posj = " << posj << '\n';
        //trovo A
        posA = massa[posi];
        //trovo N-Z
        Int_t posNmZ;
        if ((massa[posi] % 2) == ((NmenoZ[0] * NmenoZ[0]) % 2)) {posNmZ = NmenoZ[2 * posj];}
        else {posNmZ = NmenoZ[2 * posj + 1];}
        //cout << "posNmZ = " << posNmZ << '\n';
        posZ = (posA - posNmZ) / 2;
}