// Created by Chathuranga Sirimanna (chathuranga.sirimanna@wayne.edu), Gojko Vujanovic (vujanovic@wayne.edu) 
// Requires FinalState hadrons status, E, Px, Py, Pz, eta, phi as input, also requires pTHatCrossSection (in mb) as input
// Output is a txt file which contains p_T, Single hadron yield and statistical error

#include <iostream>
#include <fstream>
#include <cstdio>
#include <cmath>
#include "iomanip"
#include "string"
#include <sstream>

using namespace std;                                                                                                                  

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

    int StartTime = time(NULL);
   
    const double sigma_inel=70.0;//total inelastic p-p cross section at 5.02 TeV
 
    int NpTHatBin = 1; //Number of pTHat bins

    //These arrays are used to get the minimum and maximum values of current pTHat bin
    int pTHatMin[NpTHatBin] = {50};
    int pTHatMax[NpTHatBin] = {70};

    int NpTBin = 10; //Number of pT bins used in the Histogram
    double pTMin = 1; //Minimum pT in the Histogram in GeV
    double pTMax = 31; //Maximum pT in the Histogram in GeV
    double EtaCut =2.8; //Absolute value of the eta range
    double pTHad=0.0; //Use to store pT of each hadron before put into the histogram
    double pTBinW = (pTMax - pTMin)/NpTBin; //pT Bin width of the histogram if we use same size bins
    double HardCS[NpTHatBin]; //Store hard cross section for each pT hat bin
    double HardCSErr[NpTHatBin]; //Store error of the hard cross section for each pT hat bin
    double dNdpTCount_hard_ind[NpTHatBin][NpTBin];  //2D array used to store hadron (non recoil) count. index order: [ptHatBin] [Regular pt]
    double dNdpTCount_soft_ind[NpTHatBin][NpTBin];  //2D array used to store hadron count giving rise to source/sink terms to hydrodynamics caused by recoil hadrons. index order: [ptHatBin] [Regular pt]
    double CS_Err_hard_ind[NpTHatBin][NpTBin]; //2D array used to store hadron (non recoil) error. index order: [ptHatBin] [Regular pt]
    double CS_Err_soft_ind[NpTHatBin][NpTBin]; //2D array used to store hadron error stemming from hadrons that are source/sink terms to hydrodynamics. index order: [ptHatBin] [Regular pt]
    double TotalCS[NpTBin]; //Total cross section per pTHat bin
    double TotalCSErr[NpTBin]; //Total cross section error per pTHat bin
    int TotalEvents[NpTHatBin]; //Total events in each pT hat bin
    
    //cout << "\nStart the analysis\n\n" << endl;
    cout << "......Start Filing the Histograms.....\n" << endl;

    // For loop to open different pTHat bin files
    for (int k = 0; k<NpTHatBin; ++k)
    {
        stringstream HadronFile, HardCSFile;
        string HadronFilePath, HardCSFilePath;
        
        //path of the final state hadrons
        HadronFile << "hadrons";
        HadronFile << "_";
        HadronFile << pTHatMin[k];
        HadronFile << "_";
        HadronFile << pTHatMax[k];
        HadronFile << ".dat";
        HadronFilePath = HadronFile.str();
        
        //path of the Hard cross section file.
        HardCSFile << "cross_section";
        HardCSFile << "_";     
        HardCSFile << pTHatMin[k];
        HardCSFile << "_";
        HardCSFile << pTHatMax[k];
        HardCSFile << ".dat";
        HardCSFilePath = HardCSFile.str();
      
        //Open two files
        ifstream myfile, myfile2;
        myfile.open(HadronFilePath,ios::in);
        myfile2.open(HardCSFilePath,ios::in);

        //Some variables to store read data from file
        int SN=0,PID, stat;
        double Px, Py, Pz, E, Eta, Phi;

        string SN1, PID1, stat1, E1, Px1, Py1, Pz1, Eta1, Phi1;

        //Reset the bin arrays
        for(int j=0; j<NpTBin; j++){
            dNdpTCount_hard_ind[k][j] = 0;
            dNdpTCount_soft_ind[k][j] = 0;
        }
    
        // Read file using while loop
        TotalEvents[k] = 0;
        string MyString("#");
     
        while( myfile >> SN1 >> PID1 >> stat1 >> E1 >> Px1 >> Py1 >> Pz1 >> Eta1 >> Phi1)
        {
            //if the first string is # complete the event
            if(MyString.compare(SN1)==0) {
	   
                TotalEvents[k]++;
                continue;
            }

            //convert all the strings read to appropreate data type (double, int, etc.)
            PID=std::stoi(PID1);   // PDG Particle ID
            stat=std::stoi(stat1); // 0 if non recoil (hard induced) particle, 1 if recoil (soft induced ) particle, -1 sink/source left behind by recoil particles. 

            /* 4-momentum components */
            E=std::stod(E1);      
            Px=std::stod(Px1);     
            Py=std::stod(Py1);
            Pz=std::stod(Pz1);
            /* done */

            Eta=std::stod(Eta1); //pseudorapidity
             

            //select particles in a given eta range/particle id and bin them
            if(fabs(PID) == 211 || fabs(PID) == 111) //only pions
            {
               if(fabs(Eta)<EtaCut) //within the eta range
               {
                   pTHad = sqrt(pow(Px,2) + pow(Py,2)); //pT of the Hadron
            
                   int Hadbin = (int) ((pTHad - pTMin)/pTBinW); //calculate the corresponding pT bin
	    
                   //Bin all hadrons
                   if(Hadbin >= 0 && Hadbin < NpTBin && pTMin <= pTHad){
 	    
                       if(stat!=-1){
                           dNdpTCount_hard_ind[k][Hadbin]++; //hard-induced or non recoil particle if stat = 0
                       }
                       else if(stat==-1){
                           dNdpTCount_soft_ind[k][Hadbin]++; //if stat=-1, those are soft-induced sources/sinks to the hydrodynamical simulation left by recoil particles
                       }
                       else{
                           cout<<"Error, unknown status"<<endl; //generate error if something else
                       }
                   }
                }
            }
	 
        }//end while loop
    
        //Read and store Hard cross section and error
        myfile2 >> HardCS[k]>>HardCSErr[k];

        cout<<std::fixed<<std::setprecision(10);
        
        cout<<"Histogram infomation for pT hat bin "<<pTHatMin[k]<<"-"<<pTHatMax[k]<<" : "<<TotalEvents[k]<<endl;
        
        cout<<"\nSigma = "<<HardCS[k]<<"+/-"<<HardCSErr[k]<<endl;

        //loop over pTBins to calculate cross section error for a given pT hat bid and for the hard-induced as well as soft-induces (recoil and source/sink) partcles
        for(int j=0;j<NpTBin;j++)
        {
            if(dNdpTCount_hard_ind[k][j]> 0.0)//calculate cross section error and store it in CS_Err_hard_ind[k][j] or CS_Err_soft_ind[k][j]
            {
                CS_Err_hard_ind[k][j] = (dNdpTCount_hard_ind[k][j]*HardCS[k]/(TotalEvents[k]*sigma_inel*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut));
                CS_Err_hard_ind[k][j] *= sqrt( (1/dNdpTCount_hard_ind[k][j]) + pow(HardCSErr[k]/HardCS[k],2.0));

                cout<<"For pT Bin j = "<<j<<" \t Hard-induced BinContent = \t"<<dNdpTCount_hard_ind[k][j]<<"\t Scaled Value = "
                    <<(dNdpTCount_hard_ind[k][j]*HardCS[k])/(TotalEvents[k]*70*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut); //print out the if you want result.
            }
            else
            {
                CS_Err_hard_ind[k][j] = 0.0;
                cout<<"For pT Bin j = "<<j<<" \t Hard-induced BinContent = \t"<<dNdpTCount_hard_ind[k][j]<<"\t Scaled Value = "<<0.0;
            }

            if(dNdpTCount_soft_ind[k][j]> 0.0)//calculate cross section standard deviation and store it in CS_Err_hard_ind[k][j] or CS_Err_soft_ind[k][j]
            {
                CS_Err_soft_ind[k][j] = (dNdpTCount_soft_ind[k][j]*HardCS[k]/(TotalEvents[k]*sigma_inel*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut));
                CS_Err_soft_ind[k][j] *= sqrt( (1/dNdpTCount_soft_ind[k][j]) + pow(HardCSErr[k]/HardCS[k],2.0));

                cout<<"\t Soft-induced BinContent = \t"<<dNdpTCount_soft_ind[k][j]<<"\t Scaled Value = "
                    <<(dNdpTCount_soft_ind[k][j]*HardCS[k])/(TotalEvents[k]*70*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut)<<endl;
            }
            else
            {
                CS_Err_soft_ind[k][j] = 0.0;
                cout<<" \t Soft-induced BinContent = \t"<<dNdpTCount_soft_ind[k][j]<<"\t Scaled Value = "<<0.0<<endl;
            }
        }

        cout<<"Number of events for bin "<<pTHatMin[k]<<"-"<<pTHatMax[k]<<" : "<<TotalEvents[k]<<endl;
        cout<<".....pT hat bin completed.....\n"<<endl;
     
        myfile.close();
        myfile2.close();
    
    }//end pT hat bin loop
    cout<<"*****Filling Histograms Completed*****\n\n"<< endl;
    cout<<"-----Combining pT-hat bins in calculating the total cross section-----\n"<<endl;
    
    //arrays defined for the cross section calculation
    double CS_hard_ind[NpTBin],CS_soft_ind[NpTBin];
    double pT[NpTBin];
    //reset the defined arrays
    for(int j=0; j<NpTBin;j++)
    {
        CS_hard_ind[j] =0.0; //please use as temporary variables when calculating the TotalCS. This array is for convenience only.
        CS_soft_ind[j] =0.0; //please use as temporary variables when calculating the TotalCS. This array is for convenience only.
        TotalCS[j] =0.0;     //used for storing the final (total) cross section (summed over all pTHat bins). 
        TotalCSErr[j] = 0.0; //used for storing the uncertainty (or error) on the final (total) cross section. 
        pT[j] = 0.0; //only used for print-out purposes. 
    }
    
    cout<<"Initialized array's to zero\n"<<endl;
    
    //sum over pThat bins to calculate the resulting cross section and store in appropreate each pT bin
    for(int k=0;k<NpTHatBin;k++)
    {
        cout<<"For pT Hat Bin "<<pTHatMin[k]<< "-" <<pTHatMax[k]<<"\t CrossSection information is Below \n"<<endl;
        cout<<"pT\t\tHard induced cross section\t\tSoft induced cross section"<<endl;
        cout<<"\t\tCross Section\tVariance\t\tCross Section\tVariance"<<endl;
        for(int j=0; j<NpTBin;j++)
        {
            //cross section of non recoil (hard-induced) particles
            if(dNdpTCount_hard_ind[k][j]>0)
            {
                CS_hard_ind[j] = (dNdpTCount_hard_ind[k][j]*HardCS[k])/(TotalEvents[k]*sigma_inel*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut);
            }

            //cross section of source/sink terms originating from recoil (soft-induced) particles
            if(dNdpTCount_soft_ind[k][j]>0){
                CS_soft_ind[j] = (dNdpTCount_soft_ind[k][j]*HardCS[k])/(TotalEvents[k]*sigma_inel*2*M_PI*pTBinW*(pTMin + pTBinW * (j+0.5))*2.0*EtaCut);
            }

            //Total cross section = cross section of non-recoil (hard-induced) particles - cross section of particles giving source/sinks to the hydro
            if(CS_hard_ind[j] >= CS_soft_ind[j]){
               TotalCS[j] += CS_hard_ind[j] - CS_soft_ind[j];
            }

            //Calculate variance in the total cross section
            TotalCSErr[j] += pow( CS_Err_hard_ind[k][j], 2.0) + pow( CS_Err_soft_ind[k][j], 2.0);
            
            pT[j] = pTMin + pTBinW * (j+0.5);
            cout<<pT[j]<<"\t"<<CS_hard_ind[j]<<"\t"<<CS_Err_hard_ind[k][j]<<"\t\t"<<CS_soft_ind[j]<<"\t"<<CS_Err_soft_ind[k][j]<<endl;
        }
        cout<<"\n";
        cout<<"-----pT Hat bin completed-----\n\n"<<endl;
           
    }


    //Open output file and store final results
    ofstream CrossFile;
    CrossFile.open("HadronYield_AuAu_5020GeV_Colorless.txt",ios::out);

    cout<<"=====Final results for total crossSection=====\n"<<endl;
    cout<<"Hadron pT \t"<<"Obs. Cross Section\t"<<"Std Deviation"<<endl;
    for(int j=0;j<NpTBin;j++)
    {
        TotalCSErr[j] =sqrt(TotalCSErr[j]);
        cout<<pT[j]<<"\t"<<TotalCS[j]<<" \t     +/-"<<TotalCSErr[j]<<endl;
        CrossFile<<pT[j]<<"\t"<<TotalCS[j]<<"\t"<<TotalCSErr[j]<<endl;
    }
    CrossFile.close();
    
    cout<<"\n$$$$$ Analysis Completed $$$$$\n\n"<<endl;
    //Done

    int EndTime = time(NULL);
    int Hour = (EndTime-StartTime)/3600;
    int Minute = ((EndTime-StartTime)/60)-Hour*60;
    int Second = (EndTime-StartTime)-Hour*60*60 - Minute*60;
    cout<<"Programme run time = "<<Hour<<"::"<<Minute<<"::"<<Second<<endl;
      
    return 0;
}