ManagerNucleiCollisions.cpp

#include "ManagerNucleiCollisions.h"
#include "StringGeneration/NucleiCollision.h"

#include "TH1D.h"
#include "TH2D.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TTree.h"
#include "TMath.h"
#include "TRandom3.h"

#include <fstream>
#include <iostream>

using namespace std;

ManagerNucleiCollisions::ManagerNucleiCollisions()
{
    fOutputDirName = "outputs_ManagerNucleiCollisions";
    //    fOutputFileName = "testOutput.root";
    //    fOutputFile = 0x0;
    fDrawHistos = true;

    fFillTree = true;
}

// void ManagerNucleiCollisions::initOutputObjects()
// {
// }

// ManagerNucleiCollisions::ManagerNucleiCollisions(const ManagerNucleiCollisions& ) {}
ManagerNucleiCollisions::~ManagerNucleiCollisions() {}

void ManagerNucleiCollisions::generateEvents(NucleiCollision *fPtrNuclStruct, int nEvents)
{
    cout << "generating events..." << endl;

    // reset all hists
    // fHistSources->Reset();

    TH2D *fHist2D_impParVSnStrings = new TH2D("fHist2D_impParVSnStrings", ";impact parameter, fm;n strings", 400, 0, 20, 2001, -0.5, 2000.5);

    // ##### 10.07.2016 - save NucleiCollision info into the tree
    TTree *fNucleiCollisionsTree = new TTree("NucleiCollisionsTree", "NucleiCollisionsTree");
    const int NMaxStrings = 35000; // 35000; // 8000;

    Float_t fImpactParameter = 0;
    Float_t fNuclTreeRandomEventPlanePhi = 0;
    Float_t fNuclTreeNu = 0;
    Int_t fNuclTreeNumberOfStrings = 0;
    Int_t fMultFromAllStringsFictiveV0 = 0;
    Int_t fMultFromAllStringsFictiveMidEta = 0;
    //    Float_t fNuclTreeStringBoostAngle[NMaxStrings];
    //    Float_t fNuclTreeStringBoostMagn[NMaxStrings];
    Float_t fDistanceBetweenPartonsForString[NMaxStrings];
    Float_t fNuclTreeStringRadiusVectorAngle[NMaxStrings];
    Float_t fNuclTreeStringX[NMaxStrings];
    Float_t fNuclTreeStringY[NMaxStrings];
    short fStringOrigin[NMaxStrings];
    //    Bool_t fNuclTreeIsHardInteractionString[NMaxStrings];

    Int_t fNumberOfWoundedPartonsA = 0;
    Int_t fNumberOfWoundedPartonsB = 0;
    Int_t fNumberOfWoundedOnlyOncePartonsA = 0;
    Int_t fNumberOfWoundedOnlyOncePartonsB = 0;

    Int_t fNumberOfWoundedNucleonsA = 0;
    Int_t fNumberOfWoundedNucleonsB = 0;
    Int_t fNumberOfWoundedOnlyOnceNucleonsA = 0;
    Int_t fNumberOfWoundedOnlyOnceNucleonsB = 0;

    // May 2018: for core-corona:
    const int maxBinsCoreCorona = 23;
    int fStringsInCells[maxBinsCoreCorona][maxBinsCoreCorona];

    if (fFillTree)
    {
        fNucleiCollisionsTree->Branch("impactParameter", &fImpactParameter, "fImpactParameter/F");
        fNucleiCollisionsTree->Branch("randomEventPlanePhi", &fNuclTreeRandomEventPlanePhi, "fNuclTreeRandomEventPlanePhi/F");
        fNucleiCollisionsTree->Branch("nu", &fNuclTreeNu, "fNuclTreeNu/F");
        fNucleiCollisionsTree->Branch("numberOfStrings", &fNuclTreeNumberOfStrings, "fNuclTreeNumberOfStrings/I");

        fNucleiCollisionsTree->Branch("nWoundedPartonsA", &fNumberOfWoundedPartonsA, "fNumberOfWoundedPartonsA/I");
        fNucleiCollisionsTree->Branch("nWoundedPartonsB", &fNumberOfWoundedPartonsB, "fNumberOfWoundedPartonsB/I");
        fNucleiCollisionsTree->Branch("nWoundedOnlyOncePartonsA", &fNumberOfWoundedOnlyOncePartonsA, "fNumberOfWoundedOnlyOncePartonsA/I");
        fNucleiCollisionsTree->Branch("nWoundedOnlyOncePartonsB", &fNumberOfWoundedOnlyOncePartonsB, "fNumberOfWoundedOnlyOncePartonsB/I");

        fNucleiCollisionsTree->Branch("nWoundedNucleonsA", &fNumberOfWoundedNucleonsA, "fNumberOfWoundedNucleonsA/I");
        fNucleiCollisionsTree->Branch("nWoundedNucleonsB", &fNumberOfWoundedNucleonsB, "fNumberOfWoundedNucleonsB/I");
        fNucleiCollisionsTree->Branch("nWoundedOnlyOnceNucleonsA", &fNumberOfWoundedOnlyOnceNucleonsA, "fNumberOfWoundedOnlyOnceNucleonsA/I");
        fNucleiCollisionsTree->Branch("nWoundedOnlyOnceNucleonsB", &fNumberOfWoundedOnlyOnceNucleonsB, "fNumberOfWoundedOnlyOnceNucleonsB/I");

        if (1)
            fNucleiCollisionsTree->Branch("nStringsInCells", fStringsInCells, "fStringsInCells[23][23]/I");

        if (1)
        {
            //    fNucleiCollisionsTree->Branch("stringBoostAngle", fNuclTreeStringBoostAngle,"fNuclTreeStringBoostAngle[fNuclTreeNumberOfStrings]/F");
            //    fNucleiCollisionsTree->Branch("stringBoostMagn", fNuclTreeStringBoostMagn,"fNuclTreeStringBoostMagn[fNuclTreeNumberOfStrings]/F");
            fNucleiCollisionsTree->Branch("stringX", fNuclTreeStringX, "fNuclTreeStringX[fNuclTreeNumberOfStrings]/F");
            fNucleiCollisionsTree->Branch("stringY", fNuclTreeStringY, "fNuclTreeStringY[fNuclTreeNumberOfStrings]/F");
            fNucleiCollisionsTree->Branch("distanceBetweenPartonsForString", fDistanceBetweenPartonsForString, "fDistanceBetweenPartonsForString[fNuclTreeNumberOfStrings]/F");
            fNucleiCollisionsTree->Branch("stringRadiusVectorAngle", fNuclTreeStringRadiusVectorAngle, "fNuclTreeStringRadiusVectorAngle[fNuclTreeNumberOfStrings]/F");
            fNucleiCollisionsTree->Branch("stringOrigin", fStringOrigin, "fStringOrigin[fNuclTreeNumberOfStrings]/S");
            //    fNucleiCollisionsTree->Branch("isHardInteraction",fNuclTreeIsHardInteractionString,"fNuclTreeIsHardInteractionString[fNuclTreeNumberOfStrings]/O");
        }

        fNucleiCollisionsTree->Branch("multFromAllStringsFictiveV0", &fMultFromAllStringsFictiveV0, "fMultFromAllStringsFictiveV0/I");
        fNucleiCollisionsTree->Branch("multFromAllStringsFictiveMidEta", &fMultFromAllStringsFictiveMidEta, "fMultFromAllStringsFictiveMidEta/I");
    }
    // ##### event loop
    for (int iEvent = 0; iEvent < nEvents; iEvent++)
    {
        if (iEvent % (nEvents < 10000 ? 10 : 1000) == 0)
            cout << "generating " << (int)iEvent << "\r";
        cout.flush();

        // building event
        fPtrNuclStruct->setEventId(iEvent); // to be used in fPtrNuclStruct for output files (tmp?)
        fPtrNuclStruct->buildEvent();

        // need some event checks/triggers here?..
        if (!fPtrNuclStruct->isMBcollision())
        {
            iEvent--;
            continue;
        }

        fImpactParameter = fPtrNuclStruct->getImpactParameter();
        fNuclTreeNumberOfStrings = fPtrNuclStruct->getNstrings();

        if (fNuclTreeNumberOfStrings >= NMaxStrings)
        {
            cout << "AHTUNG: fNuclTreeNumberOfStrings >= NMaxStrings!!!" << endl;
            return;
        }

        fHist2D_impParVSnStrings->Fill(fImpactParameter, fNuclTreeNumberOfStrings);

        // ##### May 2018: core-corona:
        TH2I *fHist2DNstringsXY_thisEv_CoreCorona = fPtrNuclStruct->getHist2DNstringsXY_thisEv_CoreCorona();
        for (int i = 0; i < fHist2DNstringsXY_thisEv_CoreCorona->GetNbinsX(); i++) // x
        {
            for (int j = 0; j < fHist2DNstringsXY_thisEv_CoreCorona->GetNbinsY(); j++) // y
            {
                fStringsInCells[i][j] = fHist2DNstringsXY_thisEv_CoreCorona->GetBinContent(i + 1, j + 1);
                //                cout << fStringsInCells[i][j] << " "; //endl;
            }
            //            cout << endl;
        }

        // ##### fill tree
        if (fFillTree)
        {

            fNuclTreeRandomEventPlanePhi = fPtrNuclStruct->getRandomEventPlanePhi();
            fNuclTreeNu = fPtrNuclStruct->getNu();
            fMultFromAllStringsFictiveV0 = fPtrNuclStruct->getMultFromAllStringsFictiveV0();
            fMultFromAllStringsFictiveMidEta = fPtrNuclStruct->getMultFromAllStringsFictiveMidEta();
            // cout << fMultFromAllStringsFictive << endl;

            fNumberOfWoundedPartonsA = fPtrNuclStruct->getNumberOfWoundedPartonsA();
            fNumberOfWoundedPartonsB = fPtrNuclStruct->getNumberOfWoundedPartonsB();
            fNumberOfWoundedOnlyOncePartonsA = fPtrNuclStruct->getNumberOfWoundedOnlyOncePartonsA();
            fNumberOfWoundedOnlyOncePartonsB = fPtrNuclStruct->getNumberOfWoundedOnlyOncePartonsB();

            fNumberOfWoundedNucleonsA = fPtrNuclStruct->getNumberOfWoundedNucleonsA();
            fNumberOfWoundedNucleonsB = fPtrNuclStruct->getNumberOfWoundedNucleonsB();
            fNumberOfWoundedOnlyOnceNucleonsA = fPtrNuclStruct->getNumberOfWoundedOnlyOnceNucleonsA();
            fNumberOfWoundedOnlyOnceNucleonsB = fPtrNuclStruct->getNumberOfWoundedOnlyOnceNucleonsB();

            // ##### loop over strings
            for (int iString = 0; iString < fNuclTreeNumberOfStrings; iString++)
            {
                //            fNuclTreeStringBoostAngle[iString] = fPtrNuclStruct->getStringBoostAngle(iString);
                //            fNuclTreeStringBoostMagn[iString] = fPtrNuclStruct->getStringBoostMagn(iString);
                fNuclTreeStringRadiusVectorAngle[iString] = fPtrNuclStruct->getStringRadiusVectorAngle(iString);
                fNuclTreeStringX[iString] = fPtrNuclStruct->getStringXminusBover2(iString);
                fNuclTreeStringY[iString] = fPtrNuclStruct->getStringY(iString);
                fStringOrigin[iString] = fPtrNuclStruct->getStringOrigin(iString);
                fDistanceBetweenPartonsForString[iString] = fPtrNuclStruct->getDistanceBetweenPartonsForString(iString);
                //            fNuclTreeIsHardInteractionString[iString] = fPtrNuclStruct->isHardInteractionString(iString);
            }

            fNucleiCollisionsTree->Fill();
        }

    } // end of the event loop

    // draw the last generated event
    if (fDrawHistos)
    {
        fPtrNuclStruct->drawStatisticHists();
        //        drawStatHists();
    }
    fPtrNuclStruct->finalActions();

    // ##### additional histogram with info (cross-section, etc) - March 2018
    const int nBinLabels = 2;
    TH1D *fAdditionalInfo = new TH1D("fAdditionalInfo", "additional information from simulations;bin;value", nBinLabels, -0.5, nBinLabels - 0.5);
    TString strBinLabels[nBinLabels] =
        {
            //        "cross_seciton", // TUT OSHIBKA V NAPISANII!!!
            "cross_section",
            "mean_n_strings",
        };
    for (Int_t i = 1; i <= nBinLabels; i++)
        fAdditionalInfo->GetXaxis()->SetBinLabel(i, strBinLabels[i - 1].Data());

    fAdditionalInfo->SetBinContent(1, fPtrNuclStruct->getCrossSection());
    fAdditionalInfo->SetBinContent(2, fPtrNuclStruct->getMeanNstrings());

    //
    TCanvas *canv_impParVSnStrings = new TCanvas("canv_impParVSnStrings", "canvas impPar VS nStrings", 380, 280, 800, 600);
    fHist2D_impParVSnStrings->DrawCopy("colz");

    // write tree
    TFile *outFileNuclColTree = new TFile(Form("%s/nuclCollisionsTree_nEv%d.root", fOutputDirName.Data(), nEvents), "RECREATE");
    outFileNuclColTree->cd();
    if (fFillTree)
        fNucleiCollisionsTree->Write();
    fHist2D_impParVSnStrings->Write();
    fAdditionalInfo->Write();
    outFileNuclColTree->Close();
} // end of generateEvents()