// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.

// All rights not expressly granted are reserved.

//

// This software is distributed under the terms of the GNU General Public

// License v3 (GPL Version 3), copied verbatim in the file "COPYING".

//

// In applying this license CERN does not waive the privileges and immunities

// granted to it by virtue of its status as an Intergovernmental Organization

// or submit itself to any jurisdiction.

#include "TOFSimulation/Digitizer.h"

#include "DetectorsBase/GeometryManager.h"

#include "TOFSimulation/TOFSimParams.h"

#include "DetectorsRaw/HBFUtils.h"

#include "TCanvas.h"

#include "TFile.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TLeaf.h"

#include "TMath.h"

#include "TProfile2D.h"

#include "TRandom.h"

#include <algorithm>

#include <cassert>

using namespace o2::tof;

ClassImp(Digitizer);

// How data acquisition works in real data

/*

|<----------- 1 orbit ------------->|

------|-----------|-----------|-----------|------

^ ^ ^ ^ when triggers happen

|<--- latency ---|

|<- matching1->|

|<- matching2->|

|<- matching3->|

|<>| = overlap between two consecutive matching

Norbit = number of orbits elapsed

Nbunch = bunch in the current orbit (0:3563)

Ntdc = number of tdc counts within the matching window --> for 1/3 orbit (0:3649535)

raw time = trigger time (Norbit and Nbunch) - latency window + TDC(Ntdc)

*/

// What we implemented here (so far)

/*

|<----------- 1 orbit ------------->|

------|-----------|-----------|-----------|------

|<- matching1->|

|<- matching2->|

|<- matching3->|

|<>| = overlap between two consecutive matching windows

- OVERLAP between two consecutive windows: implemented, to be extensively checked

- NO LATENCY WINDOW (we manage it during raw encoding/decoding) then digits already corrected

NBC = Number of bunch since timeframe beginning = Norbit*3564 + Nbunch

Ntdc = here within the current BC -> (0:1023)

digit time = NBC*1024 + Ntdc

*/

void Digitizer::init()

{

// set first readout window in MC production getting

// orbitFirstSampled corresponds to the start of the concrete timeframe (it is set in O2DPG productions)

mReadoutWindowCurrent = uint64_t(o2::raw::HBFUtils::Instance().orbitFirstSampled) * Geo::NWINDOW_IN_ORBIT;

// method to initialize the parameters neede to digitize and the array of strip objects containing

// the digits belonging to a strip

initParameters();

for (Int_t i = 0; i < Geo::NSTRIPS; i++) {

for (Int_t j = 0; j < MAXWINDOWS; j++) {

mStrips[j].emplace_back(i);

if (j < MAXWINDOWS - 1) {

mMCTruthContainerNext[j] = &(mMCTruthContainer[j + 1]);

// mStripsNext[j] = &(mStrips[j + 1]);

}

}

}

}

//______________________________________________________________________

int Digitizer::process(const std::vector<HitType>* hits, std::vector<Digit>* digits)

{

const double max_hit_time = TOFSimParams::Instance().max_hit_time;

// hits array of TOF hits for a given simulated event

// digits passed from external to be filled, in continuous readout mode we will push it on mDigitsPerTimeFrame, final vector of digits

// printf("process event time = %f with %ld hits\n",mEventTime.getTimeNS(),hits->size());

uint64_t readoutwindow = getReadoutWindow(mEventTime.getTimeNS());

// determines the maximal readout window difference to a preceding RO which can still affect the current readout window

int max_readout_diff = int(max_hit_time * Geo::READOUTWINDOW_INV) + 1;

// early return based on events happening earlier than MAX_READOUT_DIFF away from current RO frame

if (readoutwindow < mReadoutWindowCurrent && mReadoutWindowCurrent - readoutwindow > max_readout_diff) {

return 0;

}

if (mContinuous && readoutwindow > mReadoutWindowCurrent) { // if we are moving in future readout windows flush previous ones (only for continuous readout mode)

digits->clear();

for (; mReadoutWindowCurrent < readoutwindow;) { // mReadoutWindowCurrent incremented in fillOutputContainer!!!!

fillOutputContainer(*digits); // fill all windows which are before (not yet stored) of the new current one

checkIfReuseFutureDigits();

} // close loop readout window

} // close if continuous

for (auto& hit : *hits) {

// TODO: put readout window counting/selection

// neglect very slow particles (low energy neutrons)

if (hit.GetTime() > max_hit_time) { // 1 mus

continue;

}

// discard hits arriving before the minimum readout window

auto hit_ro_window = getReadoutWindow(double(hit.GetTime()) + mEventTime.getTimeNS() /*+ Geo::LATENCYWINDOW*/);

// discard hits arriving too early

if (hit_ro_window < mReadoutWindowCurrent) {

continue;

}

processHit(hit, mEventTime.getTimeOffsetWrtBC() + Geo::LATENCYWINDOW);

} // end loop over hits

if (!mContinuous) { // fill output container per event

digits->clear();

fillOutputContainer(*digits);

}

return 0;

}

//______________________________________________________________________

Int_t Digitizer::processHit(const HitType& hit, Double_t event_time)

{

mNLastHit = 0;

Float_t pos[3] = {hit.GetX(), hit.GetY(), hit.GetZ()};

Float_t deltapos[3];

Int_t detInd[5];

Int_t detIndOtherPad[5];

Geo::getPadDxDyDz(pos, detInd, deltapos); // Get DetId and residuals

detIndOtherPad[0] = detInd[0], detIndOtherPad[1] = detInd[1],

detIndOtherPad[2] = detInd[2]; // same sector, plate, strip

Int_t otherraw = 0;

if (detInd[3] == 0) {

otherraw = 1;

}

Int_t iZshift = otherraw ? 1 : -1;

Int_t channel = Geo::getIndex(detInd);

Float_t charge = getCharge(hit.GetEnergyLoss());

// NOTE: FROM NOW ON THE TIME IS IN PS ... AND NOT IN NS

Double_t time = getShowerTimeSmeared((double(hit.GetTime()) + event_time) * 1E3 + 0.5 * Geo::TDCBIN, charge);

Float_t xLocal = deltapos[0];

Float_t zLocal = deltapos[2];

// extract trackID

auto trackID = hit.GetTrackID();

// PadId - Pad Identifier

// E | F --> PadId = 5 | 6

// A | B --> PadId = 1 | 2

// C | D --> PadId = 3 | 4

Int_t ndigits = 0; //Number of digits added

UInt_t istrip = channel / Geo::NPADS;

// check the fired PAD 1 (A)

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = 0;

mZLastShift[mNLastHit] = 0;

addDigit(channel, istrip, time, xLocal, zLocal, charge, 0, 0, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

// check PAD 2

detIndOtherPad[3] = otherraw;

detIndOtherPad[4] = detInd[4];

channel = Geo::getIndex(detIndOtherPad);

xLocal = deltapos[0]; // recompute local coordinates

if (otherraw) {

zLocal = deltapos[2] - Geo::ZPAD; // recompute local coordinates

} else {

zLocal = deltapos[2] + Geo::ZPAD;

}

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = 0;

mZLastShift[mNLastHit] = iZshift;

addDigit(channel, istrip, time, xLocal, zLocal, charge, 0, iZshift, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

// check PAD 3

detIndOtherPad[3] = detInd[3];

detIndOtherPad[4] = detInd[4] - 1;

if (detIndOtherPad[4] >= 0) {

channel = Geo::getIndex(detIndOtherPad);

xLocal = deltapos[0] + Geo::XPAD; // recompute local coordinates

zLocal = deltapos[2]; // recompute local coordinates

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = -1;

mZLastShift[mNLastHit] = 0;

addDigit(channel, istrip, time, xLocal, zLocal, charge, -1, 0, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

}

// check PAD 5

detIndOtherPad[3] = detInd[3];

detIndOtherPad[4] = detInd[4] + 1;

if (detIndOtherPad[4] < Geo::NPADX) {

channel = Geo::getIndex(detIndOtherPad);

xLocal = deltapos[0] - Geo::XPAD; // recompute local coordinates

zLocal = deltapos[2]; // recompute local coordinates

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = 1;

mZLastShift[mNLastHit] = 0;

addDigit(channel, istrip, time, xLocal, zLocal, charge, 1, 0, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

}

// check PAD 4

detIndOtherPad[3] = otherraw;

detIndOtherPad[4] = detInd[4] - 1;

if (detIndOtherPad[4] >= 0) {

channel = Geo::getIndex(detIndOtherPad);

xLocal = deltapos[0] + Geo::XPAD; // recompute local coordinates

if (otherraw) {

zLocal = deltapos[2] - Geo::ZPAD; // recompute local coordinates

} else {

zLocal = deltapos[2] + Geo::ZPAD;

}

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = -1;

mZLastShift[mNLastHit] = iZshift;

addDigit(channel, istrip, time, xLocal, zLocal, charge, -1, iZshift, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

}

// check PAD 6

detIndOtherPad[3] = otherraw;

detIndOtherPad[4] = detInd[4] + 1;

if (detIndOtherPad[4] < Geo::NPADX) {

channel = Geo::getIndex(detIndOtherPad);

xLocal = deltapos[0] - Geo::XPAD; // recompute local coordinates

if (otherraw) {

zLocal = deltapos[2] - Geo::ZPAD; // recompute local coordinates

} else {

zLocal = deltapos[2] + Geo::ZPAD;

}

if (isFired(xLocal, zLocal, charge)) {

ndigits++;

mXLastShift[mNLastHit] = 1;

mZLastShift[mNLastHit] = iZshift;

addDigit(channel, istrip, time, xLocal, zLocal, charge, 1, iZshift, detInd[3], trackID, hit.GetTime(), event_time * 1E3);

}

}

return ndigits;

}

//______________________________________________________________________

void Digitizer::addDigit(Int_t channel, UInt_t istrip, Double_t time, Float_t x, Float_t z, Float_t charge, Int_t iX, Int_t iZ,

Int_t padZfired, Int_t trackID, float geanttime, double t0)

{

// TOF digit requires: channel, time and time-over-threshold

if (mCalibApi->isOff(channel)) {

return;

}

time = getDigitTimeSmeared(time, x, z, charge); // add time smearing

charge *= getFractionOfCharge(x, z);

// tot tuned to reproduce 0.8% of orphans tot(=0)

Float_t totf = gRandom->Gaus(12. * Geo::NTOTBIN_PER_NS, 1.5 * Geo::NTOTBIN_PER_NS); // time-over-threshold

if (totf < 172) {

totf = 0;

}

int tot = int(totf);

Float_t xborder = Geo::XPAD * 0.5 - std::abs(x);

Float_t zborder = Geo::ZPAD * 0.5 - std::abs(z);

Float_t border = TMath::Min(xborder, zborder);

Float_t timewalkX = x * mTimeWalkeSlope;

Float_t timewalkZ = (z - (padZfired - 0.5) * Geo::ZPAD) * mTimeWalkeSlope;

if (border < 0) { // keep the effect only if hit out of pad

border *= -1;

Float_t extraTimeSmear = border * mTimeSlope;

time += gRandom->Gaus(mTimeDelay, extraTimeSmear);

} else {

border = 1 - border;

// if(border > 0) printf("deltat =%f\n",mTimeDelay*border*border*border);

// else printf("deltat=0\n");

// getchar();

if (border > 0) {

time += mTimeDelay * border * border * border;

}

}

time += TMath::Sqrt(timewalkX * timewalkX + timewalkZ * timewalkZ) - mTimeDelayCorr - mTimeWalkeSlope * 2;

// Decalibrate

float tsCorr = mCalibApi->getTimeDecalibration(channel, tot * Geo::TOTBIN_NS);

if (std::abs(tsCorr) > 200E3) { // accept correction up to 200 ns

LOG(error) << "Wrong de-calibration correction for ch = " << channel << ", tot = " << tot << " (Skip it)";

return;

}

mTimeLastHit[mNLastHit] = time;

mTotLastHit[mNLastHit] = tot;

mNLastHit++;

time -= tsCorr; // TODO: to be checked that "-" is correct, and we did not need "+" instead :-)

// let's move from time to bc, tdc

uint64_t nbc = (uint64_t)(time * Geo::BC_TIME_INPS_INV); // time elapsed in number of bunch crossing

//Digit newdigit(time, channel, (time - Geo::BC_TIME_INPS * nbc) * Geo::NTDCBIN_PER_PS, tot * Geo::NTOTBIN_PER_NS, nbc);

int tdc = int((time - Geo::BC_TIME_INPS * nbc) * Geo::NTDCBIN_PER_PS);

static long firstlongbc = long(o2::raw::HBFUtils::Instance().orbitFirstSampled) * o2::constants::lhc::LHCMaxBunches;

// add orbit and bc

nbc += mEventTime.toLong();

t0 += (mEventTime.toLong() - firstlongbc - Geo::LATENCYWINDOW_IN_BC) * Geo::BC_TIME_INPS;

// printf("orbit = %d -- bc = %d -- nbc = (%d) %d\n",mEventTime.orbit,mEventTime.bc, mEventTime.toLong(),nbc);

// printf("tdc = %d\n",tdc);

int lblCurrent = 0;

bool iscurrent = true; // if we are in the current readout window

Int_t isnext = -1;

Int_t isIfOverlap = -1;

if (mContinuous) {

isnext = nbc / Geo::BC_IN_WINDOW - mReadoutWindowCurrent;

isIfOverlap = (nbc - Geo::OVERLAP_IN_BC) / Geo::BC_IN_WINDOW - mReadoutWindowCurrent;

if (isnext == isIfOverlap) {

isIfOverlap = -1;

} else if (isnext < 0 && isIfOverlap >= 0) {

isnext = isIfOverlap;

isIfOverlap = -1;

} else if (isnext >= MAXWINDOWS && isIfOverlap < MAXWINDOWS) {

isnext = isIfOverlap;

isIfOverlap = MAXWINDOWS;

}

if (isnext < 0) {

LOG(error) << "error: isnext =" << isnext << "(current window = " << mReadoutWindowCurrent << ")"

<< " nbc = " << nbc << " -- event time = " << mEventTime.getTimeNS() << "\n";

return;

}

if (isnext < 0 || isnext >= MAXWINDOWS) {

lblCurrent = mFutureIevent.size(); // this is the size of mHeaderArray;

mFutureIevent.push_back(mEventID);

mFutureIsource.push_back(mSrcID);

mFutureItrackID.push_back(trackID);

// fill temporary digits array

insertDigitInFuture(channel, tdc, tot, nbc, lblCurrent);

return; // don't fill if doesn't match any available readout window

} else if (isIfOverlap == MAXWINDOWS) { // add in future digits but also in one of the current readout windows (beacuse of windows overlap)

lblCurrent = mFutureIevent.size();

mFutureIevent.push_back(mEventID);

mFutureIsource.push_back(mSrcID);

mFutureItrackID.push_back(trackID);

// fill temporary digits array

insertDigitInFuture(channel, tdc, tot, nbc, lblCurrent);

}

if (isnext) {

iscurrent = false;

}

}

//printf("add TOF digit c=%i n=%i\n",iscurrent,isnext);

std::vector<Strip>* strips;

o2::dataformats::MCTruthContainer<o2::tof::MCLabel>* mcTruthContainer;

if (iscurrent) {

strips = mStripsCurrent;

mcTruthContainer = mMCTruthContainerCurrent;

} else {

strips = mStripsNext[isnext - 1];

mcTruthContainer = mMCTruthContainerNext[isnext - 1];

}

fillDigitsInStrip(strips, mcTruthContainer, channel, tdc, tot, nbc, istrip, trackID, mEventID, mSrcID, geanttime, t0);

if (isIfOverlap > -1 && isIfOverlap < MAXWINDOWS) { // fill also a second readout window because of the overlap

if (!isIfOverlap) {

strips = mStripsCurrent;

mcTruthContainer = mMCTruthContainerCurrent;

} else {

strips = mStripsNext[isIfOverlap - 1];

mcTruthContainer = mMCTruthContainerNext[isIfOverlap - 1];

}

fillDigitsInStrip(strips, mcTruthContainer, channel, tdc, tot, nbc, istrip, trackID, mEventID, mSrcID, geanttime, t0);

}

}

//______________________________________________________________________

void Digitizer::fillDigitsInStrip(std::vector<Strip>* strips, o2::dataformats::MCTruthContainer<o2::tof::MCLabel>* mcTruthContainer, int channel, int tdc, int tot, uint64_t nbc, UInt_t istrip, Int_t trackID, Int_t eventID, Int_t sourceID, float geanttime, double t0)

{

int lblCurrent;

if (mcTruthContainer) {

lblCurrent = mcTruthContainer->getIndexedSize(); // this is the size of mHeaderArray;

}

Int_t lbl = (*strips)[istrip].addDigit(channel, tdc, tot, nbc, lblCurrent, 0, 0, geanttime, t0);

if (mcTruthContainer) {

if (lbl == lblCurrent) { // it means that the digit was a new one --> we have to add the info in the MC container

o2::tof::MCLabel label(trackID, eventID, sourceID, tdc);

mcTruthContainer->addElement(lbl, label);

} else {

o2::tof::MCLabel label(trackID, eventID, sourceID, tdc);

mcTruthContainer->addElementRandomAccess(lbl, label);

// sort the labels according to increasing tdc value

auto labels = mcTruthContainer->getLabels(lbl);

std::sort(labels.begin(), labels.end(),

[](o2::tof::MCLabel a, o2::tof::MCLabel b) { return a.getTDC() < b.getTDC(); });

}

}

}

//______________________________________________________________________

Double_t Digitizer::getShowerTimeSmeared(Double_t time, Float_t charge)

{

// add the smearing common to all the digits belongin to the same shower

return time + gRandom->Gaus(0, mShowerResolution);

}

//______________________________________________________________________

Double_t Digitizer::getDigitTimeSmeared(Double_t time, Float_t x, Float_t z, Float_t charge)

{

// add the smearing component which is indepedent for any digit even if belonging to the same shower (in case of

// multiple hits)

return time + gRandom->Gaus(0, mDigitResolution); // sqrt(33**2 + 50**2) ps = 60 ps

}

//______________________________________________________________________

Float_t Digitizer::getCharge(Float_t eDep)

{

// transform deposited energy in collected charge

Float_t adcMean = 50;

Float_t adcRms = 25;

return gRandom->Landau(adcMean, adcRms);

}

//______________________________________________________________________

Bool_t Digitizer::isFired(Float_t x, Float_t z, Float_t charge)

{

if (std::abs(x) > Geo::XPAD * 0.5 + 0.3) {

return kFALSE;

}

if (std::abs(z) > Geo::ZPAD * 0.5 + 0.3) {

return kFALSE;

}

Float_t effX = getEffX(x);

Float_t effZ = getEffZ(z);

Float_t efficiency = TMath::Min(effX, effZ);

if (gRandom->Rndm() > efficiency) {

return kFALSE;

}

return kTRUE;

}

//______________________________________________________________________

Float_t Digitizer::getEffX(Float_t x)

{

Float_t xborder = Geo::XPAD * 0.5 - std::abs(x);

if (xborder > 0) {

if (xborder > mBound1) {

return mEffCenter;

} else if (xborder > mBound2) {

return mEffBoundary1 + (mEffCenter - mEffBoundary1) * (xborder - mBound2) / (mBound1 - mBound2);

} else {

return mEffBoundary2 + (mEffBoundary1 - mEffBoundary2) * xborder / mBound2;

}

} else {

xborder *= -1;

if (xborder > mBound4) {

return 0;

} else if (xborder > mBound3) {

return mEffBoundary3 - mEffBoundary3 * (xborder - mBound3) / (mBound4 - mBound3);

} else {

return mEffBoundary2 + (mEffBoundary3 - mEffBoundary2) * xborder / mBound3;

}

}

return 0;

}

//______________________________________________________________________

Float_t Digitizer::getEffZ(Float_t z)

{

Float_t zborder = Geo::ZPAD * 0.5 - std::abs(z);

if (zborder > 0) {

if (zborder > mBound1) {

return mEffCenter;

} else if (zborder > mBound2) {

return mEffBoundary1 + (mEffCenter - mEffBoundary1) * (zborder - mBound2) / (mBound1 - mBound2);

} else {

return mEffBoundary2 + (mEffBoundary1 - mEffBoundary2) * zborder / mBound2;

}

} else {

zborder *= -1;

if (zborder > mBound4) {

return 0;

} else if (zborder > mBound3) {

return mEffBoundary3 - mEffBoundary3 * (zborder - mBound3) / (mBound4 - mBound3);

} else {

return mEffBoundary2 + (mEffBoundary3 - mEffBoundary2) * zborder / mBound3;

}

}

return 0;

}

//______________________________________________________________________

Float_t Digitizer::getFractionOfCharge(Float_t x, Float_t z) { return 1; }

//______________________________________________________________________

void Digitizer::setShowerSmearing()

{

mShowerResolution = 50; // smearing correlated for all digits of the same hit

if (mTOFresolution > mShowerResolution) {

mDigitResolution = TMath::Sqrt(mTOFresolution * mTOFresolution - mShowerResolution * mShowerResolution); // independent smearing for each digit

} else {

mShowerResolution = mTOFresolution;

mDigitResolution = 0;

}

}

//______________________________________________________________________

void Digitizer::initParameters()

{

// boundary references for interpolation of efficiency and resolution

mBound1 = 0.4; // distance from border when efficiency starts to decrese

mBound2 = 0.15; // second step in the fired pad

mBound3 = 0.55; // distance from border (not fired pad)

mBound4 = 0.9; // distance from border (not fired pad) when efficiency vanishes

// resolution parameters

mTOFresolution = TOFSimParams::Instance().time_resolution; // TOF global resolution in ps

setShowerSmearing();

mTimeSlope = 100; // ps/cm extra smearing if hit out of pad propto the distance from the border

mTimeDelay = 70; // time delay if hit out of pad

mTimeWalkeSlope = 40; // ps/cm

if (mMode == 0) {

mTimeSlope = 0;

mTimeDelay = 0;

mTimeWalkeSlope = 0;

}

mTimeDelayCorr = mTimeDelay / 3.5;

if (mShowerResolution > mTimeDelayCorr) {

mShowerResolution = TMath::Sqrt(mShowerResolution * mShowerResolution - mTimeDelayCorr * mTimeDelayCorr);

} else {

mShowerResolution = 0;

}

if (mShowerResolution > mTimeWalkeSlope * 0.8) {

mShowerResolution = TMath::Sqrt(mShowerResolution * mShowerResolution - mTimeWalkeSlope * mTimeWalkeSlope * 0.64);

} else {

mShowerResolution = 0;

}

// efficiency parameters

mEffCenter = TOFSimParams::Instance().eff_center; // efficiency in the center of the fired pad

mEffBoundary1 = TOFSimParams::Instance().eff_boundary1; // efficiency in mBound2

mEffBoundary2 = TOFSimParams::Instance().eff_boundary2; // efficiency in the pad border

mEffBoundary3 = TOFSimParams::Instance().eff_boundary3; // efficiency in mBound3

}

//______________________________________________________________________

void Digitizer::printParameters()

{

printf("Efficiency in the pad center = %f\n", mEffCenter);

printf("Efficiency in the pad border = %f\n", mEffBoundary2);

printf("Time resolution = %f ps (shower=%f, digit=%f)\n", mTOFresolution, mShowerResolution, mDigitResolution);

if (mTimeSlope > 0) {

printf("Degration resolution for pad with signal induced = %f ps/cm x border distance\n", mTimeSlope);

}

if (mTimeDelay > 0) {

printf("Time delay for pad with signal induced = %f ps\n", mTimeDelay);

}

if (mTimeWalkeSlope > 0) {

printf("Time walk ON = %f ps/cm\n", mTimeWalkeSlope);

}

}

//______________________________________________________________________

void Digitizer::runFullTestExample(const char* geo)

{

initParameters();

o2::tof::CalibTOFapi* api = new o2::tof::CalibTOFapi();

api->setTimeStamp(0);

api->readLHCphase();

api->readTimeSlewingParam();

setCalibApi(api);

test(geo);

}

//______________________________________________________________________

void Digitizer::test(const char* geo)

{

Int_t nhit = 1000000;

o2::base::GeometryManager::loadGeometry(geo);

o2::tof::HitType* hit = new o2::tof::HitType();

printParameters();

TH1F* h = new TH1F("hTime", "Time as from digitizer;time (ps);N", 100, -500, 500);

TH1F* h2 = new TH1F("hTot", "Tot as from digitizer;time (ns);N", 100, 0, 30);

TH1F* h3 = new TH1F("hNdig", "N_{digitis} distribution from one hit;N_{digits};N", 7, -0.5, 6.5);

TH1F* h4 = new TH1F("hTimeCorr", "Time correlation for double digits;#Deltat (ps)", 200, -1000, 1000);

TH1F* h5 = new TH1F("hTimeAv", "Time average for double digits;#Deltat (ps)", 200, -1000, 1000);

TH1F* hpad[3][3];

for (Int_t i = 0; i < 3; i++) {

for (Int_t j = 0; j < 3; j++) {

hpad[i][j] =

new TH1F(Form("hpad%i_%i", i, j), Form("Time as from digitizer, pad(%i,%i);time (ps);N", i, j), 100, -500, 500);

}

}

TProfile2D* hTimeWalk = new TProfile2D("hTimeWalk", "Time Walk;x (cm);z (cm)", 40, -1.25, 1.25, 40, -1.75, 1.75);

TH2F* hpadAll;

hpadAll = new TH2F("hpadAll", "all hits;x (cm);z (cm)", 40, -1.25, 1.25, 40, -1.75, 1.75);

TH2F* hpadHit[3][3];

TH2F* hpadEff[3][3];

for (Int_t i = 0; i < 3; i++) {

for (Int_t j = 0; j < 3; j++) {

hpadHit[i][j] = new TH2F(Form("hpadHit%i_%i", i, j), Form("pad(%i,%i) hits;x (cm);z (cm)", i - 1, 1 - j), 40,

-1.25, 1.25, 40, -1.75, 1.75);

hpadEff[i][j] = new TH2F(Form("hpadEff%i_%i", i, j), Form("pad(%i,%i) hits;x (cm);z (cm)", i - 1, 1 - j), 40,

-1.25, 1.25, 40, -1.75, 1.75);

}

}

Int_t det1[5] = {0, 0, 0, 1, 23};

Int_t det2[5] = {0, 0, 0, 0, 24};

Int_t det3[5] = {0, 0, 0, 1, 24};

Int_t det4[5] = {0, 0, 0, 0, 47};

Float_t pos[3], pos2[3], pos3[3], pos4[3];

o2::tof::Geo::getPos(det1, pos);

o2::tof::Geo::getPos(det2, pos2);

o2::tof::Geo::getPos(det3, pos3);

o2::tof::Geo::getPos(det4, pos4);

// Get strip center

pos[0] += pos2[0];

pos[1] += pos2[1];

pos[2] += pos2[2];

pos[0] *= 0.5;

pos[1] *= 0.5;

pos[2] *= 0.5;

Float_t mod;

Float_t vx[3];

Float_t vz[3];

// Get z versor

pos3[0] -= pos2[0];

pos3[1] -= pos2[1];

pos3[2] -= pos2[2];

mod = TMath::Sqrt(pos3[0] * pos3[0] + pos3[1] * pos3[1] + pos3[2] * pos3[2]);

vz[0] = pos3[0] / mod;

vz[1] = pos3[1] / mod;

vz[2] = pos3[2] / mod;

// Get x versor

pos4[0] -= pos2[0];

pos4[1] -= pos2[1];

pos4[2] -= pos2[2];

mod = TMath::Sqrt(pos4[0] * pos4[0] + pos4[1] * pos4[1] + pos4[2] * pos4[2]);

vx[0] = pos4[0] / mod;

vx[1] = pos4[1] / mod;

vx[2] = pos4[2] / mod;

Float_t x[3], dx, dz, xlocal, zlocal;

for (Int_t i = 0; i < nhit; i++) {

dx = gRandom->Rndm() * 2.5 * 48;

dz = gRandom->Rndm() * 3.5 * 2;

xlocal = dx - Int_t(dx / 2.5) * 2.5 - 1.25;

zlocal = dz - Int_t(dz / 3.5) * 3.5 - 1.75;

dx -= 2.5 * 24;

dz -= 3.5;

x[0] = pos[0] + vx[0] * dx + vz[0] * dz;

x[1] = pos[1] + vx[1] * dx + vz[1] * dz;

x[2] = pos[2] + vx[2] * dx + vz[2] * dz;

Int_t detCur[5];

o2::tof::Geo::getDetID(x, detCur);

hit->SetTime(0); // t->GetLeaf("o2root.TOF.TOFHit.mTime")->GetValue(j));

hit->SetXYZ(x[0], x[1], x[2]);

hit->SetEnergyLoss(0.0001);

processHit(*hit, mEventTime.getTimeOffsetWrtBC());

Int_t ndigits = mNLastHit;

h3->Fill(ndigits);

hpadAll->Fill(xlocal, zlocal);

for (Int_t k = 0; k < ndigits; k++) {

if (k == 0) {

h->Fill(getTimeLastHit(k));

}

if (k == 0) {

h2->Fill(getTotLastHit(k));

}

if (k == 0 && getXshift(k) == 0 && getZshift(k) == 0) {

hTimeWalk->Fill(xlocal, zlocal * (0.5 - detCur[3]) * 2, getTimeLastHit(k));

}

hpad[getXshift(k) + 1][-getZshift(k) + 1]->Fill(getTimeLastHit(k));

hpadHit[getXshift(k) + 1][-getZshift(k) + 1]->Fill(xlocal, zlocal);

}

// check double digits case (time correlations)

if (ndigits == 2) {

h4->Fill(getTimeLastHit(0) - getTimeLastHit(1));

h5->Fill((getTimeLastHit(0) + getTimeLastHit(1)) * 0.5);

}

}

h->Draw();

new TCanvas();

h2->Draw();

new TCanvas();

h3->Draw();

new TCanvas();

h4->Draw();

new TCanvas();

h5->Draw();

new TCanvas();

hTimeWalk->Draw("SURF");

TCanvas* cpad = new TCanvas();

cpad->Divide(3, 3);

for (Int_t i = 0; i < 3; i++) {

for (Int_t j = 0; j < 3; j++) {

cpad->cd(j * 3 + i + 1);

hpad[i][j]->Draw();

}

}

TCanvas* cpadH = new TCanvas();

cpadH->Divide(3, 3);

for (Int_t i = 0; i < 3; i++) {

for (Int_t j = 0; j < 3; j++) {

cpadH->cd(j * 3 + i + 1);

hpadHit[i][j]->Draw("colz");

if (j != 1) {

hpadHit[i][j]->Scale(2);

}

hpadEff[i][j]->Divide(hpadHit[i][j], hpadAll, 1, 1, "B");

hpadEff[i][j]->Draw("surf");

hpadEff[i][j]->SetMaximum(1);

hpadEff[i][j]->SetMinimum(0);

hpadEff[i][j]->SetStats(0);

}

}

printf("\nEfficiency = %f\n", (h3->GetEntries() - h3->GetBinContent(1)) / h3->GetEntries());

printf("Multiple digits fraction = %f\n\n",

(h3->GetEntries() - h3->GetBinContent(1) - h3->GetBinContent(2)) / (h3->GetEntries() - h3->GetBinContent(1)));

}

//______________________________________________________________________

void Digitizer::testFromHits(const char* geo, const char* hits)

{

o2::base::GeometryManager::loadGeometry(geo);

TFile* fHit = new TFile(hits);

fHit->ls();

TTree* t = (TTree*)fHit->Get("o2sim");

Int_t nev = t->GetEntriesFast();

o2::tof::HitType* hit = new o2::tof::HitType();

printParameters();

TH1F* h = new TH1F("hTime", "Time as from digitizer;time (ps);N", 100, -500, 500);

TH1F* h2 = new TH1F("hTot", "Tot as from digitizer;time (ns);N", 100, 0, 30);

TH1F* h3 = new TH1F("hNdig", "N_{digitis} distribution from one hit;N_{digits};N", 7, -0.5, 6.5);

for (Int_t i = 0; i < nev; i++) {

t->GetEvent(i);

Int_t nhit = t->GetLeaf("o2root.TOF.TOFHit_")->GetLen();

for (Int_t j = 0; j < nhit; j++) {

hit->SetTime(0); // t->GetLeaf("o2root.TOF.TOFHit.mTime")->GetValue(j));

hit->SetXYZ(t->GetLeaf("o2root.TOF.TOFHit.mPos.fCoordinates.fX")->GetValue(j),

t->GetLeaf("o2root.TOF.TOFHit.mPos.fCoordinates.fY")->GetValue(j),

t->GetLeaf("o2root.TOF.TOFHit.mPos.fCoordinates.fZ")->GetValue(j));

hit->SetEnergyLoss(t->GetLeaf("o2root.TOF.TOFHit.mELoss")->GetValue(j));

Int_t ndigits = processHit(*hit, mEventTime.getTimeOffsetWrtBC());

h3->Fill(ndigits);

for (Int_t k = 0; k < ndigits; k++) {

h->Fill(getTimeLastHit(k));

h2->Fill(getTotLastHit(k));

}

}

}

h->Draw();

new TCanvas();

h2->Draw();

new TCanvas();

h3->Draw();

}

//______________________________________________________________________

void Digitizer::fillOutputContainer(std::vector<Digit>& digits)

{

if (mContinuous) {

digits.clear();

mMCTruthOutputContainer->clear();

} else { // for continuos filled below

// printf("TOF fill output container\n");

// filling the digit container doing a loop on all strips

for (auto& strip : *mStripsCurrent) {

strip.fillOutputContainer(digits);

if (strip.getNumberOfDigits()) {

LOG(debug) << "strip size = " << strip.getNumberOfDigits() << " - digit size = " << digits.size() << "\n";

}

}

}

if (mContinuous) {

int first = mDigitsPerTimeFrame.size();

//printf("%i) # TOF digits = %lu (%p)\n", mIcurrentReadoutWindow, digits.size(), mStripsCurrent);

ReadoutWindowData info(first, first);

int orbit_shift = mReadoutWindowData.size() / 3;

int bc_shift = (mReadoutWindowData.size() % 3) * Geo::BC_IN_WINDOW;

info.setBCData(mFirstIR.orbit + orbit_shift, mFirstIR.bc + bc_shift);

mDigitsPerTimeFrame.insert(mDigitsPerTimeFrame.end(), digits.begin(), digits.end());

// fill diagnostics

mCalibApi->resetTRMErrors();

mCalibApi->resetDRMErrors();

float p = gRandom->Rndm();

if (mCalibApi->getEmptyTOFProb() > p) { // check empty TOF

for (int i = 0; i < Geo::kNCrate; i++) {

info.setEmptyCrate(i);

}

} else { // check empty crates when TOF is not empty

int itrmreached = -1;

bool isEmptyCrate[Geo::kNCrate];

const float* crateProb = mCalibApi->getEmptyCratesProb();

for (int i = 0; i < Geo::kNCrate; i++) {

p = gRandom->Rndm();

if (crateProb[i] > p) {

info.setEmptyCrate(i);

isEmptyCrate[i] = true;

} else { // check if filling diagnostic (noisy will be masked in clusterization, then skip here)

// Fill DRM RDH errors

for (int ie = 0; ie < mCalibApi->N_DRM_ERRORS; ie++) {

p = gRandom->Rndm();

if (mCalibApi->getDRMprobError(i, ie) > p) {

mCalibApi->processErrorDRM(i, ie);

}

}

isEmptyCrate[i] = false;

int slotreached = -1;

const std::vector<std::pair<int, float>>& trmProg = mCalibApi->getTRMerrorProb();

const std::vector<int>& trmErr = mCalibApi->getTRMmask();

for (int itrm = itrmreached + 1; itrm < trmProg.size(); itrm++) { // trm ordered by crate and slot

int crate = trmProg[itrm].first / 100;

if (crate == i) {

int slot = trmProg[itrm].first % 100;

if (slot != slotreached) { // first diagnostic of this TRM, get the random value to be compared with probability

p = gRandom->Rndm();

slotreached = slot;

}

// add diagnostic if needed

if (trmProg[itrm].second > p) {

// fill diagnostic

mCalibApi->processError(crate, slot, trmErr[itrm]);

mPatterns.push_back(slot + 28); // add slot

info.addedDiagnostic(crate);

uint32_t cbit = 1;

for (int ibit = 0; ibit < 28; ibit++) {

if (trmErr[itrm] & cbit) {

mPatterns.push_back(ibit); // add bit error

info.addedDiagnostic(crate);

}

cbit <<= 1;

}

p = 10; // no other errors allowed for this slot

} else {

p -= trmProg[itrm].second; // reduce the error probability for this slot for the next error in the same slot (sum of prob for all errors <= 1)

}

itrmreached = itrm;

} else {

break; // move to next crate

}

}

}

}

// fill strip of non-empty crates

for (auto& strip : *mStripsCurrent) {

std::map<ULong64_t, o2::tof::Digit>& dmap = strip.getDigitMap();

std::vector<ULong64_t> keyToBeRemoved;

for (auto [key, dig] : dmap) {

int crate = Geo::getCrateFromECH(Geo::getECHFromCH(dig.getChannel()));

if (isEmptyCrate[crate] || mCalibApi->isChannelError(dig.getChannel()) || mCalibApi->isChannelDRMError(dig.getChannel())) {

// flag digits to be removed

keyToBeRemoved.push_back(key);

}

}

for (auto& key : keyToBeRemoved) {

dmap.erase(key);

}

strip.fillOutputContainer(digits);

}

}

info.setNEntries(digits.size());

if (digits.size()) {

mDigitsPerTimeFrame.insert(mDigitsPerTimeFrame.end(), digits.begin(), digits.end());

}

mReadoutWindowData.push_back(info);

}

// if(! digits.size()) return;

// copying the transient labels to the output labels (stripping the tdc information)

if (mMCTruthOutputContainer) {

// copy from transientTruthContainer to mMCTruthAray

// a brute force solution for the moment; should be handled by a dedicated API

for (int index = 0; index < mMCTruthContainerCurrent->getIndexedSize(); ++index) {

mMCTruthOutputContainer->addElements(index, mMCTruthContainerCurrent->getLabels(index));

}

}

if (mContinuous) {

mMCTruthOutputContainerPerTimeFrame.push_back(*mMCTruthOutputContainer);