#include "hstarthitssimulator.h"
#include "hbeamtime.h"
#include "hruntimedb.h"
#include "hcategorymanager.h"
#include "hlinearcategory.h"
#include "hgeomvolume.h"
#include "hspecgeompar.h"
#include "hdatasource.h"

#include "hstartdef.h"
#include "hstart2cal.h"
#include "hstart2hit.h"

#include "hgeantkine.h"

#include "TFitResult.h"
#include "TRegexp.h"
#include "TH1D.h"
#include "TF1.h"
#include "TRandom.h"

#include <iostream>
#include <cstdlib>
#include <set>

//*-- Author: S.Spies
//_HADES_CLASS_DESCRIPTION
/////////////////////////////////////////////////////////////
//  HStartHitsSimulor
//
//  This class emulates StartHits to take care of delta electrons
//  from beam ions in heavy ion runs.
//
/////////////////////////////////////////////////////////////

ClassImp(HStartHitsSimulator)

//--------------------------------------------------------------------------------
    HStartHitsSimulator::HStartHitsSimulator(const Text_t* name, const Text_t* title)
    : HReconstructor(name, title),
    fCatStart2Cal(NULL), fCatStart2Hit(NULL), fCatGeantMdc(NULL), fCatGeantKine(NULL), fSpecGeomPar(NULL)
{

}

//--------------------------------------------------------------------------------
Double_t HStartHitsSimulator::FitPMF(Double_t* x, Double_t* p) {
    Int_t iMax = TMath::LocMax(fNStartHitsMax+1, p);

    for (Int_t i = 1; i <= fNStartHitsMax; ++i)
        if ((i <= iMax && p[i] < p[i-1]) || (i > iMax && p[i] > p[i-1]))
            return 0.; // Up to the maximum the probabilities must be monotonically rising and afterwards monotonically decreasing, otherwise it is no valid PMF!

    Int_t iRef = TMath::Max(0, TMath::Min(Int_t(fNStartHitsReferenceMax), TMath::FloorNint(x[0])));
    Double_t sum = 0.;
    for (Int_t i = 0; i <= fNStartHitsMax; ++i)
        sum += tmdBinomial(iRef, i) * p[i];

    return sum;
}

//--------------------------------------------------------------------------------
Bool_t HStartHitsSimulator::init() {
    if (!gHades) {
        Error("init", "No Hades object found!");
        return kFALSE;
    }

    HBeamTime::checkBeamTimeID(gHades->getBeamTimeID(),"HStartHitsSimulator::init()",kTRUE,kTRUE);

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    if (gHades->getBeamTimeID() == HADES::kApr12) {
        fStartTimeRangeReduction = (-30. - -275.)*0.04 + 45. + (80. - 15.)*0.13 + (350. - 80.)*0.13; // The reference time cut is not included as it does not influence delta electrons!
        fStartTimeReferenceMin   =  -200.;
        fStartTimeReferenceMax   =   800.;
        fStartTimeReferenceRange = fStartTimeReferenceMax - fStartTimeReferenceMin - 20. - fStartTimeRangeReduction + (-200. - -275.)*0.0125; // Here the reference time cut is included because it influences the PMF!

        if (gHades->getEmbeddingMode() > 0) { // Embedding into data
            fStartTimeMin            = -1200.;
            fStartTimeMax            =  -200.;
            fStartTimeRange          = fStartTimeMax - fStartTimeMin - (-200. - -275.)*0.0125;
            fDeltaTimeMin            = -1200.;
            fDeltaTimeMax            =  -200.;
            fDeltaTimeRange          = fDeltaTimeMax - fDeltaTimeMin - (-200. - -275.)*0.0125;
            fDeltaAddTimeMin         =  -200.;
            fDeltaAddTimeMax         =   800.;
            fDeltaAddTimeRange       = fDeltaAddTimeMax - fDeltaAddTimeMin - fStartTimeRangeReduction + (-200. - -275.)*0.0125;
        } else { // Pure Simulations
            fStartTimeMin            = -1200.;
            fStartTimeMax            =   800.;
            fStartTimeRange          = fStartTimeMax - fStartTimeMin - fStartTimeRangeReduction;
            fDeltaTimeMin            = -1200.;
            fDeltaTimeMax            =   800.;
            fDeltaTimeRange          = fDeltaTimeMax - fDeltaTimeMin - fStartTimeRangeReduction;
            fDeltaAddTimeMin         =     0.;
            fDeltaAddTimeMax         =     0.;
            fDeltaAddTimeRange       = fDeltaAddTimeMax - fDeltaAddTimeMin;
        }

        fDeltaScaleFactor        = 2.0;
        fDeltaAddScaleFactor     = 0.6;

        fNStartHitsReferenceMax  = 25;
        Double_t fNStartHitsReferencePMFTemp[fNStartHitsReferenceMax+1] = {0.2296039818, 0.1958449098, 0.1531942476, 0.1164062555, 0.0872002348, 0.0643741870, 0.0471337606, 0.0338387690, 0.0239684191, 0.0165666383, 0.0114321482, 0.0075502672, 0.0049262309, 0.0031832239, 0.0019711904, 0.0011910912, 0.0007065105, 0.0004147731, 0.0002398180, 0.0001205634, 0.0000648628, 0.0000362127, 0.0000175973, 0.0000085805, 0.0000045084, 0.0000010180}; // Day 108 gen10
        fNStartHitsReferencePMF.Set(fNStartHitsReferenceMax+1, fNStartHitsReferencePMFTemp);
        fDeltaNStartScalingBase  = 0.066;
        fDeltaNStartScalingSlope = 1.180;

        fNStartHitsMax           = 50;

        fNDeltaFiles             = 10; // Would be better to determine based on HGeantMergeSource settings, but the required functions are missing there!
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    else if (gHades->getBeamTimeID() == HADES::kMar19) {
        fStartTimeRangeReduction = (-50. - -250.)*0.03 + 70. + (350. - 20.)*0.20;
        fStartTimeReferenceMin   =  -450.;
        fStartTimeReferenceMax   =   550.;
        fStartTimeReferenceRange = fStartTimeReferenceMax - fStartTimeReferenceMin - fStartTimeRangeReduction;

        if (gHades->getEmbeddingMode() > 0) { // Embedding into data
            fStartTimeMin            = -1200.;
            fStartTimeMax            =  -450.;
            fStartTimeRange          = fStartTimeMax - fStartTimeMin;
            fDeltaTimeMin            = -1200.;
            fDeltaTimeMax            =  -200.;
            fDeltaTimeRange          = fDeltaTimeMax - fDeltaTimeMin;
            fDeltaAddTimeMin         =  -200.;
            fDeltaAddTimeMax         =   550.;
            fDeltaAddTimeRange       = fDeltaAddTimeMax - fDeltaAddTimeMin - fStartTimeRangeReduction;
        } else { // Pure Simulations
            fStartTimeMin            = -1200.;
            fStartTimeMax            =   550.;
            fStartTimeRange          = fStartTimeMax - fStartTimeMin - fStartTimeRangeReduction;
            fDeltaTimeMin            = -1200.;
            fDeltaTimeMax            =   550.;
            fDeltaTimeRange          = fDeltaTimeMax - fDeltaTimeMin - fStartTimeRangeReduction;
            fDeltaAddTimeMin         =     0.;
            fDeltaAddTimeMax         =     0.;
            fDeltaAddTimeRange       = fDeltaAddTimeMax - fDeltaAddTimeMin;
        }

        fDeltaScaleFactor        = 1.6;
        fDeltaAddScaleFactor     = 0.6;

        fNStartHitsReferenceMax  = 25;
        Double_t fNStartHitsReferencePMFTemp[fNStartHitsReferenceMax+1] = {0.2488267354, 0.2609175345, 0.1938974430, 0.1249931810, 0.0752278784, 0.0430584839, 0.0243699707, 0.0133902578, 0.0072269584, 0.0038471635, 0.0020600739, 0.0010765734, 0.0005504218, 0.0002854534, 0.0001431720, 0.0000745920, 0.0000302821, 0.0000146957, 0.0000058415, 0.0000025989, 0.0000011563, 0.0000005145, 0.0000002289, 0.0000001018, 0.0000000453, 0.0000000202}; // Day 076
        fNStartHitsReferencePMF.Set(fNStartHitsReferenceMax+1, fNStartHitsReferencePMFTemp);
        fDeltaNStartScalingBase  = 0.22;
        fDeltaNStartScalingSlope = 1.02;

        fNStartHitsMax           = 50;

        fNDeltaFiles             = 4; // Would be better to determine based on HGeantMergeSource settings, but the required functions are missing there!
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    else {
        Error("init", "Beamtime %d has no settings yet!", gHades->getBeamTimeID());
        return kFALSE;
    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    Double_t NStartClusterMean = 0.;
    for (UShort_t i = 1; i <= fNStartHitsReferenceMax; ++i)
        NStartClusterMean += i * fNStartHitsReferencePMF[i];

    fDeltaFileScaling = 1. / (NStartClusterMean / fStartTimeReferenceRange * (400. - -950.)); // 1 over average amount of Start Cluster in -950 - 400ns!

    if (gHades->getBeamTimeID() == HADES::kMar19) // Because for mar19 AgAg data AuAu delta electron simulations are used
        fDeltaFileScaling *= 47. * 47. / 79. / 79.; // Z(Ag)^2 / Z(Au)^2

    //--------------------------------------------------------------------------------
    if (gHades->getEmbeddingMode() > 0) {
        fCatStart2Cal = gHades->getCurrentEvent()->getCategory(catStart2Cal);
        if (!fCatStart2Cal)
            Warning("init", "Start2Cal category not available!");
    }

    fCatStart2Hit = gHades->getCurrentEvent()->getCategory(catStart2Hit);
    if (!fCatStart2Hit)
        Warning("init", "Start2Hit category not available!");

    fCatGeantMdc = gHades->getCurrentEvent()->getCategory(catMdcGeantRaw);
    if (!fCatGeantMdc)
        Warning("init", "GeantMdc category not available!");

    fCatGeantKine = gHades->getCurrentEvent()->getCategory(catGeantKine);
    if (!fCatGeantKine)
        Warning("init", "GeantKine category not available!");

    fSpecGeomPar = (HSpecGeomPar*) ((HRuntimeDb*) gHades->getRuntimeDb())->getContainer("SpecGeomPar");
    if (!fSpecGeomPar)
        Warning("init", "SpecGeomPar container not available!");

    //--------------------------------------------------------------------------------
    tmdBinomial    .ResizeTo(fNStartHitsReferenceMax+1, fNStartHitsMax+1);
    tmdCorrelations.ResizeTo(fNStartHitsReferenceMax+1, fNStartHitsReferenceMax+1);

    //--------------------------------------------------------------------------------
    TMatrixD tmdRefPMF(1, fNStartHitsReferenceMax+1, fNStartHitsReferencePMF.GetArray());
    Double_t p = fStartTimeReferenceRange / (((gHades->getEmbeddingMode() > 0) ? fStartTimeReferenceRange : 0.) + fStartTimeRange);
    for (Int_t iRef = 0; iRef <= fNStartHitsReferenceMax; ++iRef) {
        for (Int_t i = 0; i <= fNStartHitsMax; ++i) {
            if (iRef > i) tmdBinomial(iRef, i) = 0.; // Impossible that the smaller interval has more hits than the larger
            else          tmdBinomial(iRef, i) = TMath::Binomial(i, iRef) * TMath::Power(p, iRef) * TMath::Power(1. - p, i - iRef); // Binomial Probability
        }
    }

    TMatrixD tmdPMF = tmdRefPMF * tmdBinomial;
    tmdPMF *= 1./tmdPMF.Sum();
    fNStartHitsPMF.Set(fNStartHitsMax+1, tmdPMF.GetMatrixArray()); // Copy results

    // The result is now close but not yet perfect - A fit will do the job!
    TH1D* hRef = new TH1D("hRef", "", fNStartHitsReferenceMax+1, 0., fNStartHitsReferenceMax+1.);
    for (Int_t b = 1; b <= fNStartHitsReferenceMax+1; ++b) {
        hRef->SetBinContent(b, fNStartHitsReferencePMF[b-1]);
        hRef->SetBinError(b, TMath::Sqrt(fNStartHitsReferencePMF[b-1] * 5e6) / 5e6); // Rough estimation
    }

    TF1* fFitPMF = new TF1("fFitPMF", this, &HStartHitsSimulator::FitPMF, 0., fNStartHitsReferenceMax+1., fNStartHitsMax+1, "HStartHitsSimulator", "FitPMF");
    for (Int_t i = 0; i < fNStartHitsMax+1; ++i) {
        fFitPMF->SetParameter(i, fNStartHitsPMF[i]);
        fFitPMF->SetParError(i, fNStartHitsPMF[i] / 10.);
        fFitPMF->SetParLimits(i, 0., 1.);
    }

    TFitResultPtr tfpr = hRef->Fit(fFitPMF, "QNS");
    if (tfpr->Chi2() > 200.) { // Good Fit < 100, Bad Fit > 10000
        Error("execute", "Fit of PMF yielded large Chi2! Needs to be checked!");
        return kFALSE;
    }

    // The fit does not know that the parameters need to be normalized so I do it myself.
    Double_t norm = 0.;
    for (Int_t i = 0; i < fNStartHitsMax+1; ++i)
        norm += fFitPMF->GetParameter(i);
    for (Int_t i = 0; i < fNStartHitsMax+1; ++i)
        fNStartHitsPMF[i] = fFitPMF->GetParameter(i) / norm;

    delete hRef;
    delete fFitPMF;

    // Calculate the correlation matrix between the two intervals of interest
    if (gHades->getEmbeddingMode() > 0) {
        for (Int_t iRef = 0; iRef < fNStartHitsReferenceMax+1; ++iRef) {
            Double_t sum = 0.;

            for (Int_t i = 0; i < fNStartHitsReferenceMax+1; ++i) {
                if (iRef + i >= fNStartHitsMax+1) tmdCorrelations(iRef, i) = 0.;
                else tmdCorrelations(iRef, i) = fNStartHitsPMF[iRef+i] * TMath::Binomial(iRef+i, iRef) * TMath::Power(p, iRef) * TMath::Power(1. - p, i); // Binomial Probability
                sum += tmdCorrelations(iRef, i);
            }

            for (Int_t i = 0; i < fNStartHitsReferenceMax+1; ++i)
                tmdCorrelations(iRef, i) /= sum;
        }
    }

    return kTRUE;
}

//--------------------------------------------------------------------------------
Bool_t HStartHitsSimulator::reinit() {
    fTargetZPos.clear();

    for (Int_t t = 0; t < fSpecGeomPar->getNumTargets(); ++t) {
        HGeomVolume* tr = fSpecGeomPar->getTarget(t);
        Double_t zV1    = tr->getPoint(0)->getZ();
        Double_t zV2    = tr->getPoint(2)->getZ();
        fTargetZPos.push_back(tr->getTransform().getTransVector().getZ() + (zV1 + zV2) / 2.);
    }

    return kTRUE;
}

//--------------------------------------------------------------------------------
UShort_t HStartHitsSimulator::getRandomNStartHits(UShort_t NStartHitsReference) {
    UShort_t NStartHits = 0u;
    Double_t rndm = gRandom->Rndm();

    if (gHades->getEmbeddingMode() > 0)
        while (NStartHits <= fNStartHitsReferenceMax && tmdCorrelations(NStartHitsReference, NStartHits) <= rndm)
            rndm -= tmdCorrelations(NStartHitsReference, NStartHits++);
    else
        while (NStartHits <= fNStartHitsMax && fNStartHitsPMF[NStartHits] <= rndm)
            rndm -= fNStartHitsPMF[NStartHits++];

    return NStartHits;
}

//--------------------------------------------------------------------------------
Double_t HStartHitsSimulator::getRandomDeltaScaling() {
    if (gHades->getBeamTimeID() == HADES::kApr12) {
        if (gRandom->Rndm() < 0.16) return 0.10;
        else                        return 1.00;
    }
    else if (gHades->getBeamTimeID() == HADES::kMar19) {
        if (gRandom->Rndm() < 0.02) return 0.00;
        else                        return 1.00;
    }

    return 1.00;
}

//--------------------------------------------------------------------------------
Int_t HStartHitsSimulator::execute() {
    getAllStartHits();
    getAllStartCluster(2.);

    UShort_t NStartCluster = TMath::Min(26, getNStartCluster(2., fStartTimeReferenceMin, fStartTimeReferenceMax) - (fStartTimeReferenceMin < 0. && fStartTimeReferenceMax > 0.)); // Subtract the primary Start Cluster if the refereance range includes 0.
    UShort_t NAddStartHits = getRandomNStartHits(NStartCluster);

    if (gHades->getBeamTimeID() == HADES::kMar19 && gHades->getEmbeddingMode() > 0) { // Don't add start hits to mar19 events with extended start time window!
        TString FileName = gHades->getDataSource()->getCurrentFileName();
        if (FileName.Contains('/'))
            FileName = FileName.Data() + FileName.Last('/') + 1;

        Long64_t TimeNumber = 0;

        if (!FileName.Contains(TRegexp("^[a-z][a-z][0-9][0-9][0-3][0-9][0-9][0-5][0-9][0-5][0-9][0-5][0-9][0-9][0-9].*$")))
            Warning("execute", "Could not determine Time from FileName %s!", gHades->getDataSource()->getCurrentFileName());
        else
            TimeNumber = TString(FileName(2, 11)).Atoll();

        if (TimeNumber != 0 && TimeNumber <= 19078171630)
            NAddStartHits = 0u;
    }

    //--------------------------------------------------------------------------------
    Float_t StartTimeOffset = 0.;

    if (gHades->getEmbeddingMode() > 0) {
        if (gHades->getBeamTimeID() == HADES::kApr12)
            StartTimeOffset = gRandom->Rndm() * -40. + 10.; // Because the edge at -200ns is not sharp! (Linear rise from -230 to -190ns even though it is in Event-Time!)
        else if (gHades->getBeamTimeID() == HADES::kMar19)
            StartTimeOffset = gRandom->Rndm() * -10. - 20.; // Because the edge at -450ns is not sharp! (Linear rise from -480 to -470ns even though it is in Event-Time!)

        if (*sAllStartCluster.begin() < fStartTimeMax + StartTimeOffset)
            StartTimeOffset = *sAllStartCluster.begin() - fStartTimeMax; // Avoid overlaps
    }

    //--------------------------------------------------------------------------------
    HStart2Hit* hs2h = (HStart2Hit*) fCatStart2Hit->getObject(0);
    if (!hs2h) {
        Error("execute", "Pointer to HStart2Hit object == NULL, returning");
        return kSkipEvent;
    }

    Float_t ReactionTime = hs2h->getTime();

    //--------------------------------------------------------------------------------
    Float_t AddStartTimes[NAddStartHits];
    for (Int_t i = 0, index; i < NAddStartHits; ++i) {
        Bool_t valid = kTRUE;
        do {
            AddStartTimes[i] = gRandom->Rndm() * (fStartTimeMax - fStartTimeMin) + fStartTimeMin + StartTimeOffset; // Additional Start Times are in Event-Time!

            if (gHades->getBeamTimeID() == HADES::kApr12) {
                // Event Selection is being done in Reaction-Time!
                if (AddStartTimes[i] - ReactionTime >= -275. && AddStartTimes[i] - ReactionTime < -30.)
                    valid = gRandom->Rndm() < (1. - (AddStartTimes[i] - ReactionTime + 275.) / 245. * 0.08); // Downscaling because of Event Selection
                else if (AddStartTimes[i] - ReactionTime >= -30. && AddStartTimes[i] - ReactionTime < 15.)
                    valid = kFALSE;                                                                    // -30 +15ns basically nothing because of Event Selection
                else if (AddStartTimes[i] - ReactionTime >= 15. && AddStartTimes[i] - ReactionTime < 80.)
                    valid = gRandom->Rndm() < 0.87;                                                    // Downscaling because of Event Selection
                else if (AddStartTimes[i] - ReactionTime >= 80. && AddStartTimes[i] - ReactionTime < 350.)
                    valid = gRandom->Rndm() < (0.85 + (AddStartTimes[i] - ReactionTime - 80.) / 270. * 0.04); // Downscaling because of Event Selection
                else
                    valid = kTRUE;                                                                     // Default true
            }
            else if (gHades->getBeamTimeID() == HADES::kMar19) {
                // Event Selection is being done in Reaction-Time!
                if (AddStartTimes[i] - ReactionTime >= -250. && AddStartTimes[i] - ReactionTime < -50.)
                    valid = gRandom->Rndm() < 0.97;                                          // Downscaling because of Event Selection
                else if (AddStartTimes[i] - ReactionTime >= -50. && AddStartTimes[i] - ReactionTime < 20.)
                    valid = kFALSE;                                                          // -50 - 20ns basically nothing because of Event Selection
                else if (AddStartTimes[i] - ReactionTime >= 20. && AddStartTimes[i] - ReactionTime < 350.)
                    valid = gRandom->Rndm() < 0.80;                                          // Downscaling because of Event Selection
                else
                    valid = kTRUE;                                                           // Default true
            }

            for (Int_t j = 0; j < i && valid; ++j)
                valid = TMath::Abs(AddStartTimes[i] - AddStartTimes[j]) > 2.; // Time is in Event-Time!
        } while (!valid);

        HStart2Hit* pHit = NULL;
        if (fCatStart2Hit->getEntries() <= i+1)
            pHit = HCategoryManager::newObject(pHit, fCatStart2Hit, index);
        else
            pHit = (HStart2Hit*) fCatStart2Hit->getObject(i+1);

        pHit->setTimeAndWidth(AddStartTimes[i], 1); // Time is in Event-Time!
    }

    //--------------------------------------------------------------------------------
    std::map<Float_t, Float_t> mDelta, mDeltaAdd;
    Float_t DeltaScalingSum = 0., DeltaAddScalingSum = 0.;
    Short_t TargetLayer = getTargetLayer();

    if (TargetLayer < 0) {
        Warning("execute", "Skipped event %d because the Target Layer was not found!", gHades->getCurrentEvent()->getHeader()->getEventSeqNumber());
        return kSkipEvent;
    }

    DeltaScalingSum    += (mDelta   [-9999.] = fDeltaNStartScalingBase / 1000. * fDeltaTimeRange);    // Uncorrelated Delta Electrons
    DeltaAddScalingSum += (mDeltaAdd[-9999.] = fDeltaNStartScalingBase / 1000. * fDeltaAddTimeRange); // Uncorrelated Delta Electrons

    for (Int_t i = 0; i < NAddStartHits; ++i) {
        if (AddStartTimes[i] >= fDeltaTimeMin && AddStartTimes[i] < fDeltaTimeMax)
            DeltaScalingSum += (mDelta[AddStartTimes[i]] = getRandomDeltaScaling());
        if (AddStartTimes[i] >= fDeltaAddTimeMin && AddStartTimes[i] < fDeltaAddTimeMax)
            DeltaAddScalingSum += (mDeltaAdd[AddStartTimes[i]] = getRandomDeltaScaling());
    }

    for (std::set<Float_t>::iterator it = sAllStartCluster.begin(); it != sAllStartCluster.end(); ++it) {
        if (*it >= fDeltaTimeMin && *it < fDeltaTimeMax) {
            if (TMath::Abs(*it - ReactionTime) <= 2.) DeltaScalingSum += (mDelta[ReactionTime] = (TargetLayer + 0.5) / fTargetZPos.size());
            else                                      DeltaScalingSum += (mDelta[*it] = getRandomDeltaScaling());
        }
        if (*it >= fDeltaAddTimeMin && *it < fDeltaAddTimeMax) {
            if (TMath::Abs(*it - ReactionTime) <= 2.) DeltaAddScalingSum += (mDeltaAdd[ReactionTime] = (TargetLayer + 0.5) / fTargetZPos.size());
            else                                      DeltaAddScalingSum += (mDeltaAdd[*it] = getRandomDeltaScaling());
        }
    }

    if (DeltaScalingSum > 0.)
        for (std::map<Float_t, Float_t>::iterator it = mDelta.begin(); it != mDelta.end(); ++it)
            it->second /= DeltaScalingSum;
    if (DeltaAddScalingSum > 0.)
        for (std::map<Float_t, Float_t>::iterator it = mDeltaAdd.begin(); it != mDeltaAdd.end(); ++it)
            it->second /= DeltaAddScalingSum;

    //--------------------------------------------------------------------------------
    std::map<Int_t, std::vector<HGeantKine*> > mhgk;
    std::vector<HGeantMdc*>                    vhgm;

    for (Int_t k = 0; k < fCatGeantKine->getEntries(); ++k) {
        HGeantKine* hgk = (HGeantKine*) fCatGeantKine->getObject(k);

        if (hgk->getGeneratorInfo() == -3) { // Is a Delta Electron or a descendent of one
            if (hgk->isPrimary()) {
                if (hgk->getID() == 3)
                    mhgk[hgk->getTrack()].push_back(hgk);
            } else {
                HGeantKine* hgkPrimary = HGeantKine::getPrimary(hgk->getTrack(), (HLinearCategory*) fCatGeantKine);
                if (hgkPrimary->getID() == 3)
                    mhgk[hgkPrimary->getTrack()].push_back(hgk);
            }
        }
    }

    //--------------------------------------------------------------------------------
    if (mhgk.empty()) { // Only happens when there are no delta electrons present, so skip the event!
        Warning("execute", "Skipped event %d because it contains no delta electrons!", gHades->getCurrentEvent()->getHeader()->getEventSeqNumber());
        return kSkipEvent;
    }

    //--------------------------------------------------------------------------------
    Float_t xcoord, ycoord, tof, ptot, tOffset = -9999.;
    Double_t DeltaScaleFactor = DeltaScalingSum * fDeltaScaleFactor, DeltaAddScaleFactor = DeltaAddScalingSum * fDeltaAddScaleFactor;

    Double_t DeltaCentScaleFactor = 1.;

    for (std::map<Int_t, std::vector<HGeantKine*> >::iterator it1 = mhgk.begin(); it1 != mhgk.end(); ++it1) {
        if (gRandom->Rndm() >= (DeltaScaleFactor + DeltaAddScaleFactor) * DeltaCentScaleFactor * fDeltaNStartScalingSlope * fDeltaFileScaling / fNDeltaFiles) {
            for (std::vector<HGeantKine*>::iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
                (*it2)->setMomentum(0., 0., 0.); // Delta Electrons will be kicked out by momentum cut.
        }
        else {
            Double_t rndm = gRandom->Rndm();
            if (gRandom->Rndm() < DeltaScaleFactor / (DeltaScaleFactor + DeltaAddScaleFactor)) {
                for (std::map<Float_t, Float_t>::iterator it = mDelta.begin(); it != mDelta.end() && rndm >= 0.; ++it) {
                    tOffset = it->first - ReactionTime; // Time is in Reaction-Time! Yes, this is indeed correct because the offset is the time of the Beam Ion in Rection-Time because the simulated particles always come from 0ns!
                    rndm -= it->second;
                }
                if (tOffset < -5000.)
                    tOffset = gRandom->Rndm() * (fDeltaTimeMax - fDeltaTimeMin) + fDeltaTimeMin + StartTimeOffset - ReactionTime; // Some proportion of Delta Electrons cannot be entangled to a Start Hit
            } else {
                Double_t rndm = gRandom->Rndm();
                for (std::map<Float_t, Float_t>::iterator it = mDeltaAdd.begin(); it != mDeltaAdd.end() && rndm >= 0.; ++it) {
                    tOffset = it->first - ReactionTime; // Time is in Reaction-Time! Yes, this is indeed correct because the offset is the time of the Beam Ion in Rection-Time because the simulated particles always come from 0ns!
                    rndm -= it->second;
                }
                if (tOffset < -5000.)
                    tOffset = gRandom->Rndm() * (fDeltaAddTimeMax - fDeltaAddTimeMin) + fDeltaAddTimeMin + StartTimeOffset - ReactionTime; // Some proportion of Delta Electrons cannot be entangled to a Start Hit
            }

            for (std::vector<HGeantKine*>::iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2) {
                (*it2)->getMdcHits(vhgm);
                for (std::vector<HGeantMdc*>::iterator it3 = vhgm.begin(); it3 != vhgm.end(); ++it3) {
                    (*it3)->getHit(xcoord, ycoord, tof, ptot);
                    (*it3)->setHit(xcoord, ycoord, tof + tOffset, ptot);
                }
            }
        }
    }

    return 0;
}

//--------------------------------------------------------------------------------
void HStartHitsSimulator::getAllStartHits() {
    sAllStartHits.clear();

    if (gHades->getEmbeddingMode() == 0) return;

    HStart2Hit* hs2h = (HStart2Hit*) fCatStart2Hit->getObject(0);
    Float_t time = hs2h ? hs2h->getTime() : 0.;

    for (Int_t s = 0; s < fCatStart2Cal->getEntries(); ++s) {
        HStart2Cal* hs2c = (HStart2Cal*) fCatStart2Cal->getObject(s);

        if (hs2c->getModule() == 0 || hs2c->getModule() == 1) // Only StartX and StartY
            for (Byte_t t = 1; t <= Byte_t(TMath::Min(hs2c->getMultiplicity(), hs2c->getMaxMultiplicity())); ++t)
                if (gHades->getBeamTimeID() != HADES::kApr12 || hs2c->getTime(t) < 520. || hs2c->getTime(t) >= 540.) // Reference peak cut
                    sAllStartHits.insert(hs2c->getTime(t) - time);
    }
}

//--------------------------------------------------------------------------------
void HStartHitsSimulator::getAllStartCluster(Float_t ClusterTimeDiff) {
    sAllStartCluster.clear();

    if (sAllStartHits.empty()) return;

    UShort_t NCur = 0u;
    Float_t prevTime = -9999., SumCur = 0.;

    for (std::set<Float_t>::iterator it = sAllStartHits.begin(); it != sAllStartHits.end(); ++it) { // Iterators of a set are automatically sorted!
        if (TMath::Abs(*it - prevTime) > ClusterTimeDiff && NCur > 0u) {
            sAllStartCluster.insert(SumCur / NCur);
            SumCur = NCur = 0u;
        }

        SumCur += (prevTime = *it);
        NCur += 1u;
    }

    sAllStartCluster.insert(SumCur / NCur);
}

//--------------------------------------------------------------------------------
UShort_t HStartHitsSimulator::getNStartCluster(Float_t ClusterTimeDiff, Float_t TimeMin, Float_t TimeMax) {
    if (sAllStartHits.empty()) // Empty? No Start Hits!
        return 0u;

    UShort_t NStartCluster = 0u; // Not empty so at least one (the last in list)
    Float_t prevTime = -9999.;

    for (std::set<Float_t>::iterator it = sAllStartHits.begin(); it != sAllStartHits.end(); ++it) { // Iterators of a set are automatically sorted!
        if      (*it <  TimeMin) continue;
        else if (*it >= TimeMax) break;

        if (TMath::Abs(*it - prevTime) > ClusterTimeDiff || prevTime == -9999.)
            NStartCluster += 1u;

        prevTime = *it;
    }

    return NStartCluster;
}

//--------------------------------------------------------------------------------
Short_t HStartHitsSimulator::getTargetLayer() {
    if (fTargetZPos.empty())
        return -1;

    for (Int_t k = 0; k < fCatGeantKine->getEntries(); ++k) {
        HGeantKine* hgk = (HGeantKine*) fCatGeantKine->getObject(k);

        if (!hgk->isPrimary() || hgk->getGeneratorInfo() == -3)
            continue;

        Float_t dummy, z;
        hgk->getVertex(dummy, dummy, z);

        UShort_t closest = 0;
        for (UShort_t i = 1; i < fTargetZPos.size(); ++i)
            if (TMath::Abs(fTargetZPos[i] - z) < TMath::Abs(fTargetZPos[closest] - z))
                closest = i;

        return closest;
    }

    return -1;
}