diff --git a/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h b/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h index 05d5fdb741dbd..32e026e19d216 100644 --- a/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h +++ b/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h @@ -12,6 +12,11 @@ #ifndef DETECTORS_HMPID_BASE_INCLUDE_HMPIDDATAFORMAT_CLUSTER_H_ #define DETECTORS_HMPID_BASE_INCLUDE_HMPIDDATAFORMAT_CLUSTER_H_ +#include "CommonDataFormat/InteractionRecord.h" +#include "CommonDataFormat/RangeReference.h" +#include "DataFormatsHMP/Digit.h" +#include "HMPIDBase/Param.h" + namespace o2 { namespace hmpid @@ -21,39 +26,98 @@ namespace hmpid class Cluster { public: - Cluster() = default; - - Cluster(int chamber, int size, int NlocMax, float QRaw, float Q, float X, float Y); - ~Cluster() = default; - - int getCh() const { return mChamber; } - void setCh(int chamber) { mChamber = chamber; } - - int getSize() const { return mSize; } - void setSize(int size) { mSize = size; } + enum EClusterStatus { kFrm, + kCoG, + kLo1, + kUnf, + kMax, + kNot, + kEdg, + kSi1, + kNoLoc, + kAbn, + kBig, + kEmp = -1 }; // status flags - int getQRaw() const { return mQRaw; } - void setQRaw(int QRaw) { mQRaw = QRaw; } + public: + Cluster() : mCh(-1), mSi(-1), mSt(kEmp), mBox(-1), mNlocMax(-1), mMaxQpad(-1), mMaxQ(-1), mQRaw(0), mQ(0), mErrQ(-1), mXX(0), mErrX(-1), mYY(0), mErrY(-1), mChi2(-1), mParam(o2::hmpid::Param::instanceNoGeo()) { mDigs.clear(); }; + // Cluster(int chamber, int size, int NlocMax, float QRaw, float Q, float X, float Y) + // : mCh(chamber), mSi(size), mNlocMax(NlocMax), mQRaw(QRaw), mQ(Q), mXX(X), mYY(Y) { }; - int getQ() const { return mQ; } - void setQ(int Q) { mQ = Q; } + // Methods + // void draw(Option_t *opt=""); //overloaded TObject::Print() to draw cluster in current canvas + void print(Option_t* opt = "") const; // overloaded TObject::Print() to print cluster info + static void fitFunc(int& iNpars, double* deriv, double& chi2, double* par, int iflag); // fit function to be used by MINUIT + void coG(); // calculates center of gravity + void corrSin(); // sinoidal correction + void digAdd(o2::hmpid::Digit* pDig); // add new digit to the cluster + o2::hmpid::Digit* dig(int i) { return mDigs.at(i); } // pointer to i-th digi + inline bool isInPc(); // check if is in the current PC + void reset(); // cleans the cluster + // void setClusterParams(float xL, float yL, int iCh); //Set AliCluster3D part + int solve(std::vector* pCluLst, float* pSigmaCut, bool isUnfold); // solve cluster: MINUIT fit or CoG + // Getters + int box() { return mBox; } // Dimension of the cluster + int ch() { return mCh; } // chamber number + int size() { return mSi; } // returns number of pads in formed cluster + int status() { return mSt; } // Status of cluster + float qRaw() { return mQRaw; } // raw cluster charge in QDC channels + float q() { return mQ; } // given cluster charge in QDC channels + float qe() { return mErrQ; } // Error in cluster charge in QDC channels + float x() { return mXX; } // cluster x position in LRS + float xe() { return mErrX; } // cluster charge in QDC channels + float y() { return mYY; } // cluster y position in LRS + float ye() { return mErrY; } // cluster charge in QDC channels + float chi2() { return mChi2; } // chi2 of the fit + // Setters + void doCorrSin(bool doCorrSin) { fgDoCorrSin = doCorrSin; } // Set sinoidal correction + void setX(float x) { mXX = x; } + void setY(float y) { mYY = y; } + void setQ(float q) + { + mQ = q; + if (mQ > 4095) { + mQ = 4095; + } + } + void setQRaw(float qRaw) + { + mQRaw = qRaw; + if (mQRaw > 4095) { + mQRaw = 4095; + } + } + void setSize(int size) { mSi = size; } + void setCh(int chamber) { mCh = chamber; } + void setChi2(float chi2) { mChi2 = chi2; } + void setStatus(int status) { mSt = status; } + void findClusterSize(int i, float* pSigmaCut); // Find the clusterSize of deconvoluted clusters + virtual void clear(const Option_t*) { delete[] & mDigs; } // ===> delete [] fParam; fParam=0; } - int getX() const { return mX; } - void setX(int X) { mX = X; } + // public: + protected: + int mCh; // chamber number + int mSi; // size of the formed cluster from which this cluster deduced + int mSt; // flag to mark the quality of the cluster + int mBox; // box contaning this cluster + int mNlocMax; // number of local maxima in formed cluster + int mMaxQpad; // abs pad number of a pad with the highest charge + double mMaxQ; // that max charge value + double mQRaw; // QDC value of the raw cluster + double mQ; // QDC value of the actual cluster + double mErrQ; // error on Q + double mXX; // local x postion, [cm] + double mErrX; // error on x postion, [cm] + double mYY; // local y postion, [cm] + double mErrY; // error on y postion, [cm] + double mChi2; // some estimator of the fit quality + std::vector mDigs; //! list of digits forming this cluster - int getY() const { return mY; } - void setY(int Y) { mY = Y; } + public: + static bool fgDoCorrSin; // flag to switch on/off correction for Sinusoidal to cluster reco + Param* mParam; //! Pointer to AliHMPIDParam - protected: - int mChamber; /// chamber number - int mSize; /// size of the formed cluster from which this cluster deduced - int mNlocMax; /// number of local maxima in formed cluster - float mQRaw; /// QDC value of the raw cluster - float mQ; /// QDC value of the actual cluster - float mX; /// local x postion, [cm] - float mY; /// local y postion, [cm] - - ClassDefNV(Cluster, 1); + ClassDefNV(Cluster, 2); }; } // namespace hmpid diff --git a/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h b/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h index 926d5e9b7e1b6..b58d12ff6a259 100644 --- a/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h +++ b/DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h @@ -17,6 +17,7 @@ // History // 10/03/2021 Complete review +// 05/11/2021 Add and review for the Cluster class #ifndef DETECTORS_HMPID_BASE_INCLUDE_HMPIDDATAFORMAT_DIGIT_H_ #define DETECTORS_HMPID_BASE_INCLUDE_HMPIDDATAFORMAT_DIGIT_H_ @@ -37,11 +38,11 @@ class Digit public: // Coordinates Conversion Functions static inline uint32_t abs(int ch, int pc, int x, int y) { return ch << 24 | pc << 16 | x << 8 | y; } - static inline int ddl2C(int ddl) { return ddl >> 1; } //ddl -> chamber - static inline int a2C(uint32_t pad) { return (pad & 0xFF000000) >> 24; } //abs pad -> chamber - static inline int a2P(uint32_t pad) { return (pad & 0x00FF0000) >> 16; } //abs pad -> pc - static inline int a2X(uint32_t pad) { return (pad & 0x0000FF00) >> 8; } //abs pad -> pad X - static inline int a2Y(uint32_t pad) { return (pad & 0x000000FF); } //abs pad -> pad Y + static inline int ddl2C(int ddl) { return ddl >> 1; } // ddl -> chamber + static inline int a2C(uint32_t pad) { return (pad & 0xFF000000) >> 24; } // abs pad -> chamber + static inline int a2P(uint32_t pad) { return (pad & 0x00FF0000) >> 16; } // abs pad -> pc + static inline int a2X(uint32_t pad) { return (pad & 0x0000FF00) >> 8; } // abs pad -> pad X + static inline int a2Y(uint32_t pad) { return (pad & 0x000000FF); } // abs pad -> pad Y static inline uint32_t photo2Pad(int ch, int pc, int x, int y) { return abs(ch, pc, x, y); } static uint32_t equipment2Pad(int Equi, int Colu, int Dilo, int Chan); static uint32_t absolute2Pad(int Module, int x, int y); @@ -135,18 +136,20 @@ class Digit // pppp := photo cathode [0..5] // xxxx.xxxx := horizontal displacement [0..79] // yyyy.yyyy := vertical displacement [0..47] - //uint32_t mPad = 0; // 0xFFFFFFFF indicates invalid digit + // uint32_t mPad = 0; // 0xFFFFFFFF indicates invalid digit // Get the Geometric center of the pad - static float lorsX(int pad) { return Param::lorsX(a2P(pad), a2X(pad)); } //center of the pad x, [cm] - static float lorsY(int pad) { return Param::lorsY(a2P(pad), a2Y(pad)); } //center of the pad y, [cm] + static float lorsX(int pad) { return Param::lorsX(a2P(pad), a2X(pad)); } // center of the pad x, [cm] + static float lorsY(int pad) { return Param::lorsY(a2P(pad), a2Y(pad)); } // center of the pad y, [cm] // determines the total charge created by a hit - // might modify the localX, localY coordiates associated to the hit - static Double_t qdcTot(Double_t e, Double_t time, Int_t pc, Int_t px, Int_t py, Double_t& localX, Double_t& localY); - static Double_t intPartMathiX(Double_t x, Int_t pad); - static Double_t intPartMathiY(Double_t y, Int_t pad); - static Double_t inMathieson(Double_t localX, Double_t localY, int pad); + // might modify the localX, localY coordinates associated to the hit + static double qdcTot(double e, double time, int pc, int px, int py, double& localX, double& localY); + static double intPartMathiX(double x, int pad); + static double intPartMathiY(double y, int pad); + static double intMathieson(double localX, double localY, int pad); + static double mathiesonX(double x); // Mathieson distribution along wires X + static double mathiesonY(double x); // Mathieson distribution perp to wires Y ClassDefNV(Digit, 2); }; diff --git a/DataFormats/Detectors/HMPID/src/Cluster.cxx b/DataFormats/Detectors/HMPID/src/Cluster.cxx index 56276549e9861..cb888bb1e15ce 100644 --- a/DataFormats/Detectors/HMPID/src/Cluster.cxx +++ b/DataFormats/Detectors/HMPID/src/Cluster.cxx @@ -10,19 +10,519 @@ // or submit itself to any jurisdiction. #include +#include #include "DataFormatsHMP/Cluster.h" +#include +#include +#include ClassImp(o2::hmpid::Cluster); - namespace o2 { namespace hmpid { -Cluster::Cluster(int chamber, int size, int NlocMax, float QRaw, float Q, float X, float Y) - : mChamber(chamber), mSize(size), mNlocMax(NlocMax), mQRaw(QRaw), mQ(Q), mX(X), mY(Y) +bool o2::hmpid::Cluster::fgDoCorrSin = true; + +/*Cluster::Cluster() +{ + mCh = -1; + mSi = -1; + mSt = kEmp; + mBox = -1; + mNlocMax = -1; + mMaxQpad = -1; + mMaxQ = -1; + mQRaw = 0; + mQ = 0; + mErrQ = -1; + mXX = 0; + mErrX = -1; + mYY = 0; + mErrY = -1; + mChi2 = -1; + mDigs.clear(); + + // delete [] &mDigs; + mParam = (o2::hmpid::Param::instanceNoGeo()); + // mParam = (o2::hmpid::Param()); +}; +*/ +/*Cluster::~Cluster() +{ + delete mDigs; +} +*/ +/*void Cluster::setClusterParams(float xL,float yL,int iCh) +{ + //------------------------------------------------------------------------ + //Set the cluster properties for the AliCluster3D part + //------------------------------------------------------------------------ + mParam = o2::hmpid::Param::instance(); + if(!mParam->getInstType()) { //if there is no geometry we cannot retrieve the volId (only for monitoring) + // new(this) AliCluster3D(); return; + } + //Get the volume ID from the previously set PNEntry + uint16_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kHMPID, iCh); + + //get L->T cs matrix for a given chamber + const TGeoHMatrix *t2l = AliGeomManager::GetTracking2LocalMatrix(volId); + mParam = o2::hmpid::Param::instance(); + //transformation from the pad cs to local + xL -= 0.5 * mParam->sizeAllX(); //size of all pads with dead zones included + yL -= 0.5 * mParam->sizeAllY(); + + // Get the position in the tracking cs + float posL[3]={xL, yL, 0.}; //this is the LORS of HMPID + float posT[3]; + t2l->MasterToLocal(posL,posT); + + //Get the cluster covariance matrix in the tracking cs + float covL[9] = { 0.8 * 0.8 / 12.0, 0.0, 0.0, //pad size X + 0.0, 0.84 * 0.84 / 12.0, 0.0, //pad size Y + 0.0, 0.0, 0.1 }; //just 1 , no Z dimension ??? + + TGeoHMatrix m; + m.SetRotation(covL); + m.Multiply(t2l); + const TGeoHMatrix& t2li = t2l->Inverse(); + m.MultiplyLeft(&t2li); + float *covT = m.GetRotationMatrix(); + + // ===> new(this) AliCluster3D(volId,posT[0],posT[1],posT[2],covT[0],covT[1],covT[2],covT[4],covT[5],covT[8],0x0); // No MC labels ? +} +*/ +// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Calculates naive cluster position as a center of gravity of its digits. +void Cluster::coG() +{ + int minPadX = 999; + int minPadY = 999; + int maxPadX = -1; + int maxPadY = -1; // for box finding + + if (mDigs.size() == 0) { + return; + } // no digits in this cluster + mXX = mYY = mQRaw = 0; // init summable parameters + mCh = -1; // init chamber + int maxQpad = -1; + int maxQ = -1; // to calculate the pad with the highest charge + + o2::hmpid::Digit* pDig = nullptr; + for (int iDig = 0; iDig < mDigs.size(); iDig++) { // digits loop + int x, y, mod; + int padId = mDigs[iDig]->getPadID(); + o2::hmpid::Digit::pad2Absolute(padId, &mod, &x, &y); + if (x > maxPadX) { + maxPadX = x; + } // find the minimum box that contain the cluster MaxX + if (y > maxPadY) { + maxPadY = y; + } // MaxY + if (x < minPadX) { + minPadX = x; + } // MinX + if (y < minPadY) { + minPadY = y; + } // MinY + + float q = mDigs[iDig]->mQ; // get QDC + mXX += o2::hmpid::Digit::lorsX(padId) * q; + mYY += o2::hmpid::Digit::lorsY(padId) * q; // add digit center weighted by QDC + mQRaw += q; // increment total charge + if (q > maxQ) { + maxQpad = padId; + maxQ = (int)q; + } // to find pad with highest charge + mCh = mod; // initialize chamber number + } // digits loop + + mBox = (maxPadX - minPadX + 1) * 100 + maxPadY - minPadY + 1; // dimension of the box: format Xdim*100+Ydim + if (mQRaw != 0) { + mXX /= mQRaw; + mYY /= mQRaw; + } // final center of gravity + + if (mDigs.size() > 1 && fgDoCorrSin) { + corrSin(); + } // correct it by sinoid + + mQ = mQRaw; // Before starting fit procedure, Q and QRaw must be equal + mMaxQpad = maxQpad; + mMaxQ = maxQ; // store max charge pad to the field + mChi2 = 0; // no Chi2 to find + mNlocMax = 0; // proper status from this method + mSt = o2::hmpid::Cluster::kCoG; + return; + // setClusterParams(mXX, mYY, mCh); //need to fill the AliCluster3D part +} // CoG() + +// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Correction of cluster x position due to sinoid, see HMPID TDR page 30 +void Cluster::corrSin() +{ + int pc; + int px; + int py; + o2::hmpid::Param::lors2Pad(mXX, mYY, pc, px, py); // tmp digit to get it center + double x = mXX - mParam->lorsX(pc, px); // diff between cluster x and center of the pad contaning this cluster + mXX += 3.31267e-2 * sin(2 * M_PI / 0.8 * x) - 2.66575e-3 * sin(4 * M_PI / 0.8 * x) + 2.80553e-3 * sin(6 * M_PI / 0.8 * x) + 0.0070; + return; +} + +// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +/*void Cluster::draw(Option_t*) +{ + TMarker *pMark = new TMarker(X(), Y(), 5); + pMark->SetUniqueID(mSt);pMark->SetMarkerColor(kBlue); pMark->Draw(); +} +*/ +// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Cluster::fitFunc(int& iNpars, double* deriv, double& chi2, double* par, int iflag) +{ + // Cluster fit function + // par[0]=x par[1]=y par[2]=q for the first Mathieson shape + // par[3]=x par[4]=y par[5]=q for the second Mathieson shape and so on up to iNpars/3 Mathieson shapes + // For each pad of the cluster calculates the difference between actual pad charge and the charge induced to this pad by all Mathieson distributions + // Then the chi2 is calculated as the sum of this value squared for all pad in the cluster. + // Arguments: iNpars - number of parameters which is number of local maxima of cluster * 3 + // chi2 - function result to be minimised + // par - parameters array of size iNpars + // Returns: none + Cluster* pClu = (Cluster*)TVirtualFitter::GetFitter()->GetObjectFit(); + int nPads = pClu->mSi; + chi2 = 0; + int iNshape = iNpars / 3; + for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster + double dQpadMath = 0; + for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad + double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID()); + dQpadMath += par[3 * j + 2] * fracMathi; // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson + } + if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) { + chi2 += std::pow((pClu->dig(i)->mQ - dQpadMath), 2.0) / pClu->dig(i)->mQ; // chi2 function to be minimized + } + } + //---calculate gradients... + if (iflag == 2) { + float** derivPart; + derivPart = new float*[iNpars]; + for (int j = 0; j < iNpars; j++) { + deriv[j] = 0; + derivPart[j] = new float[nPads]; + for (int i = 0; i < nPads; i++) { + derivPart[j][i] = 0; + } + } + for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster + for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad + double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID()); + double lx = o2::hmpid::Digit::lorsX(pClu->dig(i)->getPadID()); + double ly = o2::hmpid::Digit::lorsY(pClu->dig(i)->getPadID()); + derivPart[3 * j][i] += par[3 * j + 2] * (o2::hmpid::Digit::mathiesonX(par[3 * j] - lx - 0.5 * o2::hmpid::Param::sizePadX()) - o2::hmpid::Digit::mathiesonX(par[3 * j] - lx + 0.5 * o2::hmpid::Param::sizePadX())) * + o2::hmpid::Digit::intPartMathiY(par[3 * j + 1], pClu->dig(i)->getPadID()); + derivPart[3 * j + 1][i] += par[3 * j + 2] * (o2::hmpid::Digit::mathiesonY(par[3 * j + 1] - ly - 0.5 * o2::hmpid::Param::sizePadY()) - o2::hmpid::Digit::mathiesonY(par[3 * j + 1] - ly + 0.5 * o2::hmpid::Param::sizePadY())) * + o2::hmpid::Digit::intPartMathiX(par[3 * j], pClu->dig(i)->getPadID()); + derivPart[3 * j + 2][i] += fracMathi; + } + } + // loop on all pads of the cluster + for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster + double dQpadMath = 0; // pad charge collector + for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad + float fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID()); + dQpadMath += par[3 * j + 2] * fracMathi; + if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) { + deriv[3 * j] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j][i]; + deriv[3 * j + 1] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 1][i]; + deriv[3 * j + 2] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 2][i]; + } + } + } + // delete array... + for (int i = 0; i < iNpars; i++) { + delete[] derivPart[i]; + } + delete[] derivPart; + } //---gradient calculations ended + // fit ended. Final calculations +} // FitFunction() + +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Print current cluster +void Cluster::print(Option_t* opt) const +{ + const char* status = nullptr; + switch (mSt) { + case kFrm: + status = "formed "; + break; + case kUnf: + status = "unfolded (fit)"; + break; + case kCoG: + status = "coged "; + break; + case kLo1: + status = "locmax 1 (fit)"; + break; + case kMax: + status = "exceeded (cog)"; + break; + case kNot: + status = "not done (cog)"; + break; + case kEmp: + status = "empty "; + break; + case kEdg: + status = "edge (fit)"; + break; + case kSi1: + status = "size 1 (cog)"; + break; + case kNoLoc: + status = "no LocMax(fit)"; + break; + case kAbn: + status = "Abnormal fit "; + break; + case kBig: + status = "Big Clu(>100) "; + break; + default: + status = "??????"; + break; + } + float ratio = 0; + if (mQ > 0.0 && mQRaw > 0.0) { + ratio = mQ / mQRaw * 100; + } + printf("%sCLU: ch=%i (X = %7.3f, Y = %7.3f) Q=%8.3f Qraw=%8.3f(%3.0f%%) Size=%2i DimBox=%i LocMax=%i Chi2=%7.3f %s\n", + opt, mCh, mXX, mYY, mQ, mQRaw, ratio, mSi, mBox, mNlocMax, mChi2, status); + if (mDigs.size() > 0) { + std::cout << "Digits of Cluster" << std::endl; + for (int i; i < mDigs.size(); i++) { + std::cout << mDigs.at(i) << std::endl; + } + } + return; +} // Print() + +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, bool isTryUnfold) +{ + // This methode is invoked when the cluster is formed to solve it. Solve the cluster means to try to unfold the cluster + // into the local maxima number of clusters. This methode is invoked by AliHMPIDRconstructor::Dig2Clu() on cluster by cluster basis. + // At this point, cluster contains a list of digits, cluster charge and size is precalculated in AddDigit(), position is preset to (-1,-1) in ctor, + // status is preset to kFormed in AddDigit(), chamber-sector info is preseted to actual values in AddDigit() + // Method first finds number of local maxima and if it's more then one tries to unfold this cluster into local maxima number of clusters + // Arguments: pCluLst - cluster list pointer where to add new cluster(s) + // isTryUnfold - flag to switch on/off unfolding + // Returns: number of local maxima of original cluster + const int kMaxLocMax = 6; // max allowed number of loc max for fitting + coG(); // First calculate CoG for the given cluster + int iCluCnt = pCluLst->size(); // get current number of clusters already stored in the list by previous operations + int rawSize = mSi; // get current raw cluster size + if (rawSize > 100) { + mSt = kBig; + } else if (isTryUnfold == false) { + mSt = kNot; + } else if (rawSize == 1) { + mSt = kSi1; + } + if (rawSize > 100 || isTryUnfold == false || rawSize == 1) { // No deconv if: 1 - big cluster (also avoid no zero suppression!) + // setClusterParams(mXX, mYY, mCh); // 2 - flag is set to FALSE + // new ((*pCluLst)[iCluCnt++]) Cluster(*this); // 3 - size = 1 + pCluLst->push_back(o2::hmpid::Cluster(*this)); + return 1; // add this raw cluster + } + + // Phase 0. Initialise Fitter + double arglist[10]; + float ierflg = 0.; + TVirtualFitter* fitter = TVirtualFitter::Fitter((TObject*)this, 3 * 6); // initialize Fitter + arglist[0] = -1; + ierflg = fitter->ExecuteCommand("SET PRI", arglist, 1); // no printout + ierflg = fitter->ExecuteCommand("SET NOW", arglist, 0); // no warning messages + arglist[0] = 1; + ierflg = fitter->ExecuteCommand("SET GRA", arglist, 1); // force Fitter to use my gradient + fitter->SetFCN(Cluster::fitFunc); + // Phase 1. Find number of local maxima. Strategy is to check if the current pad has QDC more then all neigbours. Also find the box contaning the cluster + mNlocMax = 0; + for (int iDig1 = 0; iDig1 < rawSize; iDig1++) { // first digits loop + o2::hmpid::Digit* pDig1 = mDigs.at(iDig1); // take next digit + int iCnt = 0; // counts how many neighbouring pads has QDC more then current one + for (int iDig2 = 0; iDig2 < rawSize; iDig2++) { // loop on all digits again + if (iDig1 == iDig2) { + continue; + } // the same digit, no need to compare + o2::hmpid::Digit* pDig2 = mDigs.at(iDig2); // take second digit to compare with the first one + int dist = TMath::Sign(int(pDig1->mX - pDig2->mX), 1) + TMath::Sign(int(pDig1->mY - pDig2->mY), 1); // distance between pads + if (dist == 1) { // means dig2 is a neighbour of dig1 + if (pDig2->mQ >= pDig1->mQ) { + iCnt++; // count number of pads with Q more then Q of current pad + } + } + } // second digits loop + if (iCnt == 0 && mNlocMax < kMaxLocMax) { // this pad has Q more then any neighbour so it's local maximum + float xStart = o2::hmpid::Digit::lorsX(pDig1->getPadID()); + float yStart = o2::hmpid::Digit::lorsY(pDig1->getPadID()); + float xMin = xStart - mParam->sizePadX(); + float xMax = xStart + mParam->sizePadX(); + float yMin = yStart - mParam->sizePadY(); + float yMax = yStart + mParam->sizePadY(); + ierflg = fitter->SetParameter(3 * mNlocMax, Form("x%i", mNlocMax), xStart, 0.1, xMin, xMax); // X,Y,Q initial values of the loc max pad + ierflg = fitter->SetParameter(3 * mNlocMax + 1, Form("y%i", mNlocMax), yStart, 0.1, yMin, yMax); // X, Y constrained to be near the loc max + ierflg = fitter->SetParameter(3 * mNlocMax + 2, Form("q%i", mNlocMax), pDig1->mQ, 0.1, 0, 10000); // Q constrained to be positive + mNlocMax++; + } // if this pad is local maximum + } // first digits loop + + // Phase 2. Fit loc max number of Mathiesons or add this current cluster to the list + // case 1 -> no loc max found + if (mNlocMax == 0) { // case of no local maxima found: pads with same charge... + mNlocMax = 1; + mSt = kNoLoc; + // setClusterParams(mXX, mYY, mCh); //need to fill the AliCluster3D part + pCluLst->push_back(o2::hmpid::Cluster(*this)); // add new unfolded cluster pCluLst->push_back(o2::hmpid::Cluster(*this)); + return mNlocMax; + } + + // case 2 -> loc max found. Check # of loc maxima + if (mNlocMax >= kMaxLocMax) { + // setClusterParams(mXX, mYY, mCh); // if # of local maxima exceeds kMaxLocMax... + mSt = kMax; + pCluLst->push_back(o2::hmpid::Cluster(*this)); //...add this raw cluster + } else { // or resonable number of local maxima to fit and user requested it + // Now ready for minimization step + arglist[0] = 500; // number of steps and sigma on pads charges + arglist[1] = 1.; // + ierflg = fitter->ExecuteCommand("SIMPLEX", arglist, 2); // start fitting with Simplex + if (!ierflg) { + fitter->ExecuteCommand("MIGRAD", arglist, 2); // fitting improved by Migrad + } + if (ierflg) { + double strategy = 2.; + ierflg = fitter->ExecuteCommand("SET STR", &strategy, 1); // change level of strategy + if (!ierflg) { + ierflg = fitter->ExecuteCommand("SIMPLEX", arglist, 2); // start fitting with Simplex + if (!ierflg) { + fitter->ExecuteCommand("MIGRAD", arglist, 2); // fitting improved by Migrad + } + } + } + if (ierflg) { + mSt = kAbn; // no convergence of the fit... + } + double dummy; + char sName[80]; // vars to get results from Minuit + double edm; + double errdef; + int nvpar; + int nparx; + for (int i = 0; i < mNlocMax; i++) { // store the local maxima parameters + fitter->GetParameter(3 * i, sName, mXX, mErrX, dummy, dummy); // X + fitter->GetParameter(3 * i + 1, sName, mYY, mErrY, dummy, dummy); // Y + fitter->GetParameter(3 * i + 2, sName, mQ, mErrQ, dummy, dummy); // Q + fitter->GetStats(mChi2, edm, errdef, nvpar, nparx); // get fit infos + if (mNlocMax > 1) { + findClusterSize(i, pSigmaCut); // find clustersize for deconvoluted clusters + // after this call, fSi temporarly is the calculated size. Later is set again + // to its original value + } + if (mSt != kAbn) { + if (mNlocMax != 1) { + mSt = kUnf; // if unfolded + } + if (mNlocMax == 1 && mSt != kNoLoc) { + mSt = kLo1; // if only 1 loc max + } + if (!isInPc()) { + mSt = kEdg; // if Out of Pc + } + if (mSt == kNoLoc) { + mNlocMax = 0; // if with no loc max (pads with same charge..) + } + } + // setClusterParams(mXX, mYY, mCh); //need to fill the AliCluster3D part + pCluLst->push_back(o2::hmpid::Cluster(*this)); // add new unfolded cluster + if (mNlocMax > 1) { + setSize(rawSize); // Original raw size is set again to its proper value + } + } + } + return mNlocMax; +} // Solve() - {}; +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// Estimate of the clustersize for a deconvoluted cluster +void Cluster::findClusterSize(int i, float* pSigmaCut) +{ + int size = 0; + for (int iDig = 0; iDig < mSi; iDig++) { // digits loop + o2::hmpid::Digit* pDig = dig(iDig); // take digit + int iCh = pDig->mCh; + double qPad = mQ * o2::hmpid::Digit::intMathieson(x(), y(), pDig->getPadID()); // pad charge pDig-> + // AliDebug(1,Form("Chamber %i X %i Y %i SigmaCut %i pad %i qpadMath %8.2f qPadRaw %8.2f Qtotal %8.2f cluster n.%i", + // iCh, o2::hmpid::Digit::a2X(pDig->getPadID()), o2::hmpid::Digit::a2Y(pDig->getPadID()), + // pSigmaCut[iCh],iDig,qPad,pDig->mQ,mQRaw,i)); + if (qPad > pSigmaCut[iCh]) { + size++; + } + } + // AliDebug(1,Form(" Calculated size %i",size)); + if (size > 0) { + setSize(size); // in case of size == 0, original raw clustersize used + } +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Bool_t Cluster::isInPc() +{ + // Check if (X,Y) position is inside the PC limits + // Arguments: + // Returns: True or False + int pc = ((o2::hmpid::Digit*)&mDigs.at(0))->getPh(); // (o2::hmpid::Digit*)&mDigs.at(iDig) + + if (mXX < Param::minPcX(pc) || mXX > Param::maxPcX(pc) || mYY < Param::minPcY(pc) || mYY > Param::maxPcY(pc)) { + return false; + } + + return true; +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Cluster::digAdd(Digit* pDig) +{ + // Adds a given digit to the list of digits belonging to this cluster, cluster is not owner of digits + // Arguments: pDig - pointer to digit to be added + // Returns: none + if (mDigs.size() == 0) { // create list of digits in the first invocation + mSi = 0; + // std::vector fDigs; + } + // fDigs->Add(pDig); + mDigs.push_back(pDig); + mSt = kFrm; + mSi++; + return; +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Cluster::reset() +{ + // + // + // + if (mDigs.size() > 0) { + mDigs.clear(); + } + // mDigs={0x0}; + mSt = kEmp; + mQRaw = mQ = 0; + mXX = mYY = 0; + mCh = mSi = mBox = mNlocMax = mMaxQpad = -1; + mMaxQ = mErrQ = mErrX = mErrY = mChi2 = -1; // empty ctor +} } // namespace hmpid } // namespace o2 diff --git a/DataFormats/Detectors/HMPID/src/Digit.cxx b/DataFormats/Detectors/HMPID/src/Digit.cxx index a654ea611ab7b..397d1a3cd5d0d 100644 --- a/DataFormats/Detectors/HMPID/src/Digit.cxx +++ b/DataFormats/Detectors/HMPID/src/Digit.cxx @@ -9,7 +9,6 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// /// \file Digit.cxx /// \author Antonio Franco - INFN Bari /// \brief Base Class to manage HMPID Digit data @@ -18,7 +17,7 @@ /* ------ HISTORY --------- 10/03/2021 / complete review - + 05/11/2021 Add and review for the Cluster class */ #include @@ -271,7 +270,7 @@ float Digit::getFractionalContributionForPad(HitType const& hit, int somepad) // converting chamber id and hit coordiates to local coordinates Param::instance()->mars2Lors(chamber, tmp, localX, localY); // calculate charge fraction in given pad - return Digit::inMathieson(localX, localY, somepad); + return Digit::intMathieson(localX, localY, somepad); } /// QdcTot : Samples total charge associated to a hit @@ -297,19 +296,19 @@ Double_t Digit::qdcTot(Double_t e, Double_t time, Int_t pc, Int_t px, Int_t py, return 0; } else { double y = Param::lorsY(pc, py); - localY = ((y - localY) > 0) ? y - 0.2 : y + 0.2; //shift to the nearest anod wire + localY = ((y - localY) > 0) ? y - 0.2 : y + 0.2; // shift to the nearest anod wire - double x = (localX > 66.6) ? localX - 66.6 : localX; //sagita is for PC (0-64) and not for chamber - double qdcEle = 34.06311 + 0.2337070 * x + 5.807476e-3 * x * x - 2.956471e-04 * x * x * x + 2.310001e-06 * x * x * x * x; //reparametrised from DiMauro + double x = (localX > 66.6) ? localX - 66.6 : localX; // sagita is for PC (0-64) and not for chamber + double qdcEle = 34.06311 + 0.2337070 * x + 5.807476e-3 * x * x - 2.956471e-04 * x * x * x + 2.310001e-06 * x * x * x * x; // reparametrised from DiMauro int iNele = int((e / 26e-9) * 0.8); if (iNele < 1) { - iNele = 1; //number of electrons created by hit, if photon e=0 implies iNele=1 + iNele = 1; // number of electrons created by hit, if photon e=0 implies iNele=1 } for (Int_t i = 1; i <= iNele; i++) { double rnd = gRandom->Rndm(); if (rnd == 0) { - rnd = 1e-12; //1e-12 is a protection against 0 from rndm + rnd = 1e-12; // 1e-12 is a protection against 0 from rndm } Q -= qdcEle * TMath::Log(rnd); } @@ -324,15 +323,15 @@ Double_t Digit::qdcTot(Double_t e, Double_t time, Int_t pc, Int_t px, Int_t py, /// @param[in] x : position of the center of Mathieson distribution /// @param[in] pad : the Digit Unique Id [0x00CPXXYY] /// @return : a charge fraction [0-1] imposed into the pad -Double_t Digit::intPartMathiX(Double_t x, int pad) +double Digit::intPartMathiX(double x, int pad) { - Double_t shift1 = -lorsX(pad) + 0.5 * o2::hmpid::Param::sizePadX(); - Double_t shift2 = -lorsX(pad) - 0.5 * o2::hmpid::Param::sizePadX(); + double shift1 = -lorsX(pad) + 0.5 * o2::hmpid::Param::sizePadX(); + double shift2 = -lorsX(pad) - 0.5 * o2::hmpid::Param::sizePadX(); - Double_t ux1 = o2::hmpid::Param::sqrtK3x() * TMath::TanH(o2::hmpid::Param::k2x() * (x + shift1) / o2::hmpid::Param::pitchAnodeCathode()); - Double_t ux2 = o2::hmpid::Param::sqrtK3x() * TMath::TanH(o2::hmpid::Param::k2x() * (x + shift2) / o2::hmpid::Param::pitchAnodeCathode()); + double ux1 = o2::hmpid::Param::sqrtK3x() * tanh(o2::hmpid::Param::k2x() * (x + shift1) / o2::hmpid::Param::pitchAnodeCathode()); + double ux2 = o2::hmpid::Param::sqrtK3x() * tanh(o2::hmpid::Param::k2x() * (x + shift2) / o2::hmpid::Param::pitchAnodeCathode()); - return o2::hmpid::Param::k4x() * (TMath::ATan(ux2) - TMath::ATan(ux1)); + return o2::hmpid::Param::k4x() * (atan(ux2) - atan(ux1)); } /// IntPartMathiY : Integration of Mathieson. @@ -342,18 +341,18 @@ Double_t Digit::intPartMathiX(Double_t x, int pad) /// @param[in] y : position of the center of Mathieson distribution /// @param[in] pad : the Digit Unique Id [0x00CPXXYY] /// @return : a charge fraction [0-1] imposed into the pad -Double_t Digit::intPartMathiY(Double_t y, int pad) +double Digit::intPartMathiY(double y, int pad) { - Double_t shift1 = -lorsY(pad) + 0.5 * o2::hmpid::Param::sizePadY(); - Double_t shift2 = -lorsY(pad) - 0.5 * o2::hmpid::Param::sizePadY(); + double shift1 = -lorsY(pad) + 0.5 * o2::hmpid::Param::sizePadY(); + double shift2 = -lorsY(pad) - 0.5 * o2::hmpid::Param::sizePadY(); - Double_t uy1 = o2::hmpid::Param::sqrtK3y() * TMath::TanH(o2::hmpid::Param::k2y() * (y + shift1) / o2::hmpid::Param::pitchAnodeCathode()); - Double_t uy2 = o2::hmpid::Param::sqrtK3y() * TMath::TanH(o2::hmpid::Param::k2y() * (y + shift2) / o2::hmpid::Param::pitchAnodeCathode()); + double uy1 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift1) / o2::hmpid::Param::pitchAnodeCathode()); + double uy2 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift2) / o2::hmpid::Param::pitchAnodeCathode()); - return o2::hmpid::Param::k4y() * (TMath::ATan(uy2) - TMath::ATan(uy1)); + return o2::hmpid::Param::k4y() * (atan(uy2) - atan(uy1)); } -/// InMathieson : Integration of Mathieson. +/// IntMathieson : Integration of Mathieson. /// This is the answer to electrostatic problem of charge distrubution in MWPC /// described elsewhere. (NIM A370(1988)602-603) /// @@ -361,10 +360,39 @@ Double_t Digit::intPartMathiY(Double_t y, int pad) /// @param[in] localY : Y position of the center of Mathieson distribution /// @param[in] pad : the Digit Unique Id [0x00CPXXYY] /// @return : a charge fraction [0-1] imposed into the pad -Double_t Digit::inMathieson(Double_t localX, Double_t localY, Int_t pad) +double Digit::intMathieson(double localX, double localY, int pad) { return 4. * intPartMathiX(localX, pad) * intPartMathiY(localY, pad); } + +// Mathieson function. +// This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603) +// Arguments: x- position of the center of Mathieson distribution +// Returns: value of the Mathieson function +double Digit::mathiesonX(double x) +{ + double lambda = x / o2::hmpid::Param::pitchAnodeCathode(); + double tanh_v = tanh(o2::hmpid::Param::k2x() * lambda); + double a = 1 - tanh_v * tanh_v; + double b = 1 + o2::hmpid::Param::sqrtK3x() * o2::hmpid::Param::sqrtK3x() * tanh_v * tanh_v; + double mathi = o2::hmpid::Param::k1x() * a / b; + return mathi; +} + +// Mathieson function. +// This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603) +// Arguments: x- position of the center of Mathieson distribution +// Returns: value of the Mathieson function +double Digit::mathiesonY(double y) +{ + double lambda = y / o2::hmpid::Param::pitchAnodeCathode(); + double tanh_v = tanh(o2::hmpid::Param::k2y() * lambda); + double a = 1 - tanh_v * tanh_v; + double b = 1 + o2::hmpid::Param::sqrtK3y() * o2::hmpid::Param::sqrtK3y() * tanh_v * tanh_v; + double mathi = o2::hmpid::Param::k1y() * a / b; + return mathi; +} + /* // ---- Time conversion functions ---- diff --git a/Detectors/HMPID/base/include/HMPIDBase/Param.h b/Detectors/HMPID/base/include/HMPIDBase/Param.h index 7466ee9f26101..9d5b044e155cf 100644 --- a/Detectors/HMPID/base/include/HMPIDBase/Param.h +++ b/Detectors/HMPID/base/include/HMPIDBase/Param.h @@ -30,7 +30,7 @@ class Param { public: - //ctor&dtor + // ctor&dtor virtual ~Param() { if (fgInstance) { @@ -42,33 +42,33 @@ class Param } } - void print(Option_t* opt = "") const; //print current parametrization + void print(Option_t* opt = "") const; // print current parametrization - static Param* instance(); //pointer to Param singleton - static inline Param* instanceNoGeo(); //pointer to Param singleton without geometry.root for MOOD, displays, ... - //geo info + static Param* instance(); // pointer to Param singleton + static Param* instanceNoGeo(); // pointer to Param singleton without geometry.root for MOOD, displays, ... + // geo info enum EChamberData { kMinCh = 0, kMaxCh = 6, kMinPc = 0, - kMaxPc = 5 }; //Segmenation + kMaxPc = 5 }; // Segmenation enum EPadxData { kPadPcX = 80, kMinPx = 0, kMaxPx = 79, - kMaxPcx = 159 }; //Segmentation structure along x + kMaxPcx = 159 }; // Segmentation structure along x enum EPadyData { kPadPcY = 48, kMinPy = 0, kMaxPy = 47, - kMaxPcy = 143 }; //Segmentation structure along y - //The electronics takes the 32bit int as: first 9 bits for the pedestal and the second 9 bits for threshold - // - values below should be within range + kMaxPcy = 143 }; // Segmentation structure along y + // The electronics takes the 32bit int as: first 9 bits for the pedestal and the second 9 bits for threshold + // - values below should be within range enum EPedestalData { kPadMeanZeroCharge = 400, kPadSigmaZeroCharge = 20, kPadMeanMasked = 401, - kPadSigmaMasked = 20 }; //One can go up to 5 sigma cut, overflow is protected in AliHMPIDCalib + kPadSigmaMasked = 20 }; // One can go up to 5 sigma cut, overflow is protected in AliHMPIDCalib static float r2d() { return 57.2957795; } - static float sizePadX() { return fgCellX; } //pad size x, [cm] - static float sizePadY() { return fgCellY; } //pad size y, [cm] + static float sizePadX() { return fgCellX; } // pad size x, [cm] + static float sizePadY() { return fgCellY; } // pad size y, [cm] static float sizePcX() { return fgPcX; } // PC size x static float sizePcY() { return fgPcY; } // PC size y @@ -76,32 +76,32 @@ class Param static float maxPcY(Int_t iPc) { return fgkMaxPcY[iPc]; } // PC limits static float minPcX(Int_t iPc) { return fgkMinPcX[iPc]; } // PC limits static float minPcY(Int_t iPc) { return fgkMinPcY[iPc]; } // PC limits - static Int_t nsig() { return fgNSigmas; } //Getter n. sigmas for noise - static float sizeAllX() { return fgAllX; } //all PCs size x, [cm] - static float sizeAllY() { return fgAllY; } //all PCs size y, [cm] + static Int_t nsig() { return fgNSigmas; } // Getter n. sigmas for noise + static float sizeAllX() { return fgAllX; } // all PCs size x, [cm] + static float sizeAllY() { return fgAllY; } // all PCs size y, [cm] - //center of the pad x, [cm] + // center of the pad x, [cm] static float lorsX(Int_t pc, Int_t padx) { return (padx + 0.5) * sizePadX() + fgkMinPcX[pc]; } - //center of the pad y, [cm] + // center of the pad y, [cm] static float lorsY(Int_t pc, Int_t pady) { return (pady + 0.5) * sizePadY() + fgkMinPcY[pc]; } - //PhiMin (degree) of the camber ch + // PhiMin (degree) of the camber ch float chPhiMin(Int_t ch) { return lors2Mars(ch, lorsX(ch, kMinPx) - mX, lorsY(ch, kMinPy) - mY).Phi() * r2d(); } - //ThMin (degree) of the camber ch + // ThMin (degree) of the camber ch float chThMin(Int_t ch) { return lors2Mars(ch, lorsX(ch, kMinPx) - mX, lorsY(ch, kMinPy) - mY).Theta() * r2d(); } - //PhiMax (degree) of the camber ch + // PhiMax (degree) of the camber ch float chPhiMax(Int_t ch) { return lors2Mars(ch, lorsX(ch, kMaxPcx) - mX, lorsY(ch, kMaxPcy) - mY).Phi() * r2d(); } - //ThMax (degree) of the camber ch + // ThMax (degree) of the camber ch float chThMax(Int_t ch) { return lors2Mars(ch, lorsX(ch, kMaxPcx) - mX, lorsY(ch, kMaxPcy) - mY).Theta() * r2d(); } static void lors2Pad(float x, float y, Int_t& pc, Int_t& px, Int_t& py); //(x,y)->(pc,px,py) - static bool isOverTh(float q) { return q >= fgThreshold; } //is digit over threshold? + static bool isOverTh(float q) { return q >= fgThreshold; } // is digit over threshold? - bool getInstType() const { return fgInstanceType; } //return if the instance is from geom or ideal + bool getInstType() const { return fgInstanceType; } // return if the instance is from geom or ideal - inline static bool isInDead(float x, float y); //is the point in dead area? - inline static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch); //is a dead pad? + inline static bool isInDead(float x, float y); // is the point in dead area? + inline static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch); // is a dead pad? inline void setChStatus(Int_t ch, bool status = kTRUE); inline void setSectStatus(Int_t ch, Int_t sect, bool status); @@ -109,7 +109,7 @@ class Param inline void printChStatus(Int_t ch); inline void setGeomAccept(); - static Int_t inHVSector(float y); //find HV sector + static Int_t inHVSector(float y); // find HV sector static Int_t radiator(float y) { if (inHVSector(y) < 0) { @@ -126,15 +126,15 @@ class Param } return y - radiator(y) * fgkMinPcY[radiator(y)]; } - //is point inside chamber boundaries? + // is point inside chamber boundaries? static bool isInside(float x, float y, float d = 0) { return x > -d && y > -d && x < fgkMaxPcX[kMaxPc] + d && y < fgkMaxPcY[kMaxPc] + d; } - //For optical properties + // For optical properties static double ePhotMin() { return 5.5; } // - static double ePhotMax() { return 8.5; } //Photon energy range,[eV] + static double ePhotMax() { return 8.5; } // Photon energy range,[eV] static double nIdxRad(double eV, double temp) { return TMath::Sqrt(1 + 0.554 * (1239.84 / eV) * (1239.84 / eV) / ((1239.84 / eV) * (1239.84 / eV) - 5769)) - 0.0005 * (temp - 20); @@ -143,16 +143,16 @@ class Param static double nMgF2Idx(double eV) { return 1.7744 - 2.866e-3 * (1239.842609 / eV) + 5.5564e-6 * (1239.842609 / eV) * (1239.842609 / eV); } // MgF2 idx of trasparency system static double nIdxGap(double eV) { return 1 + 0.12489e-6 / (2.62e-4 - eV * eV / 1239.84 / 1239.84); } static double lAbsRad(double eV) { return (eV < 7.8) * (gausPar(eV, 3.20491e16, -0.00917890, 0.742402) + gausPar(eV, 3035.37, 4.81171, 0.626309)) + (eV >= 7.8) * 0.0001; } - static double lAbsWin(double eV) { return (eV < 8.2) * (818.8638 - 301.0436 * eV + 36.89642 * eV * eV - 1.507555 * eV * eV * eV) + (eV >= 8.2) * 0.0001; } //fit from DiMauro data 28.10.03 + static double lAbsWin(double eV) { return (eV < 8.2) * (818.8638 - 301.0436 * eV + 36.89642 * eV * eV - 1.507555 * eV * eV * eV) + (eV >= 8.2) * 0.0001; } // fit from DiMauro data 28.10.03 static double lAbsGap(double eV) { return (eV < 7.75) * 6512.399 + (eV >= 7.75) * 3.90743e-2 / (-1.655279e-1 + 6.307392e-2 * eV - 8.011441e-3 * eV * eV + 3.392126e-4 * eV * eV * eV); } - static double qEffCSI(double eV) { return (eV > 6.07267) * 0.344811 * (1 - exp(-1.29730 * (eV - 6.07267))); } //fit from DiMauro data 28.10.03 + static double qEffCSI(double eV) { return (eV > 6.07267) * 0.344811 * (1 - exp(-1.29730 * (eV - 6.07267))); } // fit from DiMauro data 28.10.03 static double gausPar(double x, double a1, double a2, double a3) { return a1 * TMath::Exp(-0.5 * ((x - a2) / a3) * ((x - a2) / a3)); } - //find the temperature of the C6F14 in a given point with coord. y (in x is uniform) + // find the temperature of the C6F14 in a given point with coord. y (in x is uniform) inline static double findTemp(double tLow, double tUp, double y); double getEPhotMean() const { return mPhotEMean; } - double getRefIdx() const { return mRefIdx; } //running refractive index + double getRefIdx() const { return mRefIdx; } // running refractive index double meanIdxRad() const { return nIdxRad(mPhotEMean, mTemp); } double meanIdxWin() const { return nIdxWin(mPhotEMean); } @@ -167,10 +167,10 @@ class Param double winIdx() const { return 1.5787; } //<--TEMPORAR--> to be removed in future. Mean refractive index of WIN material (SiO2) double gapIdx() const { return 1.0005; } //<--TEMPORAR--> to be removed in future. Mean refractive index of GAP material (CH4) - static Int_t stack(Int_t evt = -1, Int_t tid = -1); //Print stack info for event and tid - static Int_t stackCount(Int_t pid, Int_t evt); //Counts stack particles of given sort in given event - static void idealPosition(Int_t iCh, TGeoHMatrix* m); //ideal position of given chamber - //trasformation methodes + static Int_t stack(Int_t evt = -1, Int_t tid = -1); // Print stack info for event and tid + static Int_t stackCount(Int_t pid, Int_t evt); // Counts stack particles of given sort in given event + static void idealPosition(Int_t iCh, TGeoHMatrix* m); // ideal position of given chamber + // trasformation methodes void lors2Mars(Int_t c, double x, double y, double* m, Int_t pl = kPc) const { double z = 0; @@ -193,14 +193,14 @@ class Param double m[3]; lors2Mars(c, x, y, m, pl); return TVector3(m); - } //MRS->LRS + } // MRS->LRS void mars2Lors(Int_t c, double* m, double& x, double& y) const { double l[3]; mM[c]->MasterToLocal(m, l); x = l[0] + mX; y = l[1] + mY; - } //MRS->LRS + } // MRS->LRS void mars2LorsVec(Int_t c, double* m, double& th, double& ph) const { double l[3]; @@ -209,37 +209,37 @@ class Param th = TMath::ATan(pt / l[2]); ph = TMath::ATan2(l[1], l[0]); } - void lors2MarsVec(Int_t c, double* m, double* l) const { mM[c]->LocalToMasterVect(m, l); } //LRS->MRS + void lors2MarsVec(Int_t c, double* m, double* l) const { mM[c]->LocalToMasterVect(m, l); } // LRS->MRS TVector3 norm(Int_t c) const { double n[3]; norm(c, n); return TVector3(n); - } //norm + } // norm void norm(Int_t c, double* n) const { double l[3] = {0, 0, 1}; mM[c]->LocalToMasterVect(l, n); - } //norm - void point(Int_t c, double* p, Int_t plane) const { lors2Mars(c, 0, 0, p, plane); } //point of given chamber plane + } // norm + void point(Int_t c, double* p, Int_t plane) const { lors2Mars(c, 0, 0, p, plane); } // point of given chamber plane - void setTemp(double temp) { mTemp = temp; } //set actual temperature of the C6F14 - void setEPhotMean(double ePhotMean) { mPhotEMean = ePhotMean; } //set mean photon energy + void setTemp(double temp) { mTemp = temp; } // set actual temperature of the C6F14 + void setEPhotMean(double ePhotMean) { mPhotEMean = ePhotMean; } // set mean photon energy - void setRefIdx(double refRadIdx) { mRefIdx = refRadIdx; } //set running refractive index + void setRefIdx(double refRadIdx) { mRefIdx = refRadIdx; } // set running refractive index - void setNSigmas(Int_t sigmas) { fgNSigmas = sigmas; } //set sigma cut - void setThreshold(Int_t thres) { fgThreshold = thres; } //set sigma cut - void setInstanceType(bool inst) { fgInstanceType = inst; } //kTRUE if from geomatry kFALSE if from ideal geometry - //For PID - double sigLoc(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); //error due to cathode segmetation - double sigGeom(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); //error due to unknown photon origin - double sigCrom(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); //error due to unknonw photon energy - double sigma2(double trkTheta, double trkPhi, double ckovTh, double ckovPh); //photon candidate sigma^2 + void setNSigmas(Int_t sigmas) { fgNSigmas = sigmas; } // set sigma cut + void setThreshold(Int_t thres) { fgThreshold = thres; } // set sigma cut + void setInstanceType(bool inst) { fgInstanceType = inst; } // kTRUE if from geomatry kFALSE if from ideal geometry + // For PID + double sigLoc(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); // error due to cathode segmetation + double sigGeom(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); // error due to unknown photon origin + double sigCrom(double trkTheta, double trkPhi, double ckovTh, double ckovPh, double beta); // error due to unknonw photon energy + double sigma2(double trkTheta, double trkPhi, double ckovTh, double ckovPh); // photon candidate sigma^2 - static double sigmaCorrFact(Int_t iPart, double occupancy); //correction factor for theoretical resolution + static double sigmaCorrFact(Int_t iPart, double occupancy); // correction factor for theoretical resolution - //Mathieson Getters + // Mathieson Getters static double pitchAnodeCathode() { return fgkD; } static double sqrtK3x() { return fgkSqrtK3x; } @@ -253,18 +253,18 @@ class Param // enum EPlaneId { kPc, kRad, - kAnod }; //3 planes in chamber + kAnod }; // 3 planes in chamber enum ETrackingFlags { kMipDistCut = -9, kMipQdcCut = -5, - kNoPhotAccept = -11 }; //flags for Reconstruction + kNoPhotAccept = -11 }; // flags for Reconstruction protected: - static /*const*/ float fgkMinPcX[6]; //limits PC - static /*const*/ float fgkMinPcY[6]; //limits PC - static /*const*/ float fgkMaxPcX[6]; //limits PC + static /*const*/ float fgkMinPcX[6]; // limits PC + static /*const*/ float fgkMinPcY[6]; // limits PC + static /*const*/ float fgkMaxPcX[6]; // limits PC static /*const*/ float fgkMaxPcY[6]; - static bool fgMapPad[160][144][7]; //map of pads to evaluate if they are active or dead (160,144) pads for 7 chambers + static bool fgMapPad[160][144][7]; // map of pads to evaluate if they are active or dead (160,144) pads for 7 chambers // Mathieson constants // For HMPID --> x direction means parallel to the wires: K3 = 0.66 (NIM A270 (1988) 602-603) fig.1 @@ -277,24 +277,24 @@ class Param static const double fgkSqrtK3y, fgkK2y, fgkK1y, fgkK4y; // - static Int_t fgNSigmas; //sigma Cut - static Int_t fgThreshold; //sigma Cut - static bool fgInstanceType; //kTRUE if from geomatry kFALSE if from ideal geometry + static Int_t fgNSigmas; // sigma Cut + static Int_t fgThreshold; // sigma Cut + static bool fgInstanceType; // kTRUE if from geomatry kFALSE if from ideal geometry - static float fgCellX, fgCellY, fgPcX, fgPcY, fgAllX, fgAllY; //definition of HMPID geometric parameters - Param(bool noGeo); //default ctor is protected to enforce it to be singleton + static float fgCellX, fgCellY, fgPcX, fgPcY, fgAllX, fgAllY; // definition of HMPID geometric parameters + Param(bool noGeo); // default ctor is protected to enforce it to be singleton - static Param* fgInstance; //static pointer to instance of Param singleton + static Param* fgInstance; // static pointer to instance of Param singleton - TGeoHMatrix* mM[7]; //pointers to matrices defining HMPID chambers rotations-translations - float mX; //x shift of LORS with respect to rotated MARS - float mY; //y shift of LORS with respect to rotated MARS - double mRefIdx; //running refractive index of C6F14 - double mPhotEMean; //mean energy of photon - double mTemp; //actual temparature of C6F14 + TGeoHMatrix* mM[7]; // pointers to matrices defining HMPID chambers rotations-translations + float mX; // x shift of LORS with respect to rotated MARS + float mY; // y shift of LORS with respect to rotated MARS + double mRefIdx; // running refractive index of C6F14 + double mPhotEMean; // mean energy of photon + double mTemp; // actual temparature of C6F14 private: - Param(const Param& r); //dummy copy constructor - Param& operator=(const Param& r); //dummy assignment operator + Param(const Param& r); // dummy copy constructor + Param& operator=(const Param& r); // dummy assignment operator ClassDefNV(Param, 1); }; diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h index 15d738a3363bf..c53ae8cab0da2 100644 --- a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h @@ -16,6 +16,7 @@ #include #include +#include "HMPIDBase/Param.h" #include "DataFormatsHMP/Cluster.h" #include "DataFormatsHMP/Digit.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -39,20 +40,27 @@ class Clusterer Clusterer(const Clusterer&) = delete; Clusterer& operator=(const Clusterer&) = delete; - void process(std::vector const& digits, std::vector& clusters, MCLabelContainer const* digitMCTruth); + // void process(std::vector const& digits, std::vector& clusters, MCLabelContainer const* digitMCTruth); - void setMCTruthContainer(o2::dataformats::MCTruthContainer* truth) { mClsLabels = truth; } + // void setMCTruthContainer(o2::dataformats::MCTruthContainer* truth) { mClsLabels = truth; } + + static void Dig2Clu(std::vector* mDigs, std::vector* mClus, float* pUserCut, bool isUnfold = kTRUE); // digits->clusters + static void FormClu(Cluster* pClu, int pDig, std::vector* pDigList, TMatrixF* pDigMap); // cluster formation recursive algorithm + static int UseDig(int padX, int padY, TMatrixF* pDigMap); // use this pad's digit to form a cluster + inline bool IsDigSurvive(Digit* pDig) const; // check for sigma cut private: - //void processChamber(std::vector& clusters, MCLabelContainer const* digitMCTruth); - //void fetchMCLabels(const Digit* dig, std::array& labels, int& nfilled) const; + // void processChamber(std::vector& clusters, MCLabelContainer const* digitMCTruth); + // void fetchMCLabels(const Digit* dig, std::array& labels, int& nfilled) const; o2::dataformats::MCTruthContainer* mClsLabels = nullptr; // Cluster MC labels - Digit* mContributingDigit[6]; //! array of digits contributing to the cluster; this will not be stored, it is temporary to build the final cluster - int mNumberOfContributingDigits; //! number of digits contributing to the cluster; this will not be stored, it is temporary to build the final cluster - void addContributingDigit(Digit* dig); - void buildCluster(Cluster& c, MCLabelContainer const* digitMCTruth); + // Digit* mContributingDigit[6]; //! array of digits contributing to the cluster; this will not be stored, it is temporary to build the final cluster + // int mNumberOfContributingDigits; //! number of digits contributing to the cluster; this will not be stored, it is temporary to build the final cluster + // std::vector mDigs; + // std::vector mClus; + // void addContributingDigit(Digit* dig); + // void buildCluster(Cluster& c, MCLabelContainer const* digitMCTruth); }; } // namespace hmpid diff --git a/Detectors/HMPID/reconstruction/src/Clusterer.cxx b/Detectors/HMPID/reconstruction/src/Clusterer.cxx index eac7c3da8f518..26b9bf7c066b0 100644 --- a/Detectors/HMPID/reconstruction/src/Clusterer.cxx +++ b/Detectors/HMPID/reconstruction/src/Clusterer.cxx @@ -14,6 +14,7 @@ #include #include "FairLogger.h" // for LOG #include "Framework/Logger.h" +#include "HMPIDBase/Param.h" #include "DataFormatsHMP/Cluster.h" #include "HMPIDReconstruction/Clusterer.h" #include "SimulationDataFormat/MCCompLabel.h" @@ -22,8 +23,134 @@ using namespace o2::hmpid; -//__________________________________________________ -void Clusterer::process(std::vector const& digits, std::vector& clusters, MCLabelContainer const* digitMCTruth) +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +/*Clusterer::Clusterer() +{ + // mDigs = {nullptr}; + // mClus = {nullptr}; +}*/ +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Clusterer::Dig2Clu(std::vector* mDigs, std::vector* mClus, float* pUserCut, bool isTryUnfold) +{ + // Finds all clusters for a given digits list provided not empty. Currently digits list is a list of all digits for a single chamber. + // Puts all found clusters in separate lists, one per clusters. + // Arguments: pDigAll - list of digits for all chambers + // pCluAll - list of clusters for all chambers + // isTryUnfold - flag to choose between CoG and Mathieson fitting + // Returns: none + + TMatrixF padMap(Param::kMinPx, Param::kMaxPcx, Param::kMinPy, Param::kMaxPcy); // pads map for single chamber 0..159 x 0..143 + + int pUsedDig = -1; + int padChX = 0, padChY = 0, module = 0; + o2::hmpid::Cluster* clu; // tmp cluster to be used as current + + for (int iCh = Param::kMinCh; iCh <= Param::kMaxCh; iCh++) { // chambers loop + padMap = (Float_t)-1; // reset map to -1 (means no digit for this pad) + for (int iDig = 0; iDig < mDigs->size(); iDig++) { // digits loop for current chamber + o2::hmpid::Digit::pad2Absolute(mDigs->at(iDig).getPadID(), &module, &padChX, &padChY); + if (module != iCh) { + continue; + } + padMap(padChX, padChY) = iDig; // fill the map for the given chamber, (padx,pady) cell takes digit index + } // digits loop for current chamber + + for (int iDig = 0; iDig < mDigs->size(); iDig++) { // digits loop for current chamber + o2::hmpid::Digit::pad2Absolute(mDigs->at(iDig).getPadID(), &module, &padChX, &padChY); + if (module != iCh) { + continue; + } + pUsedDig = UseDig(padChX, padChY, &padMap); + if (pUsedDig == -1) { + continue; + } // this digit is already taken in FormClu(), go after next digit + clu = new o2::hmpid::Cluster; + clu->setCh(iCh); + FormClu(clu, pUsedDig, mDigs, &padMap); // form cluster starting from this digit by recursion + clu->solve(mClus, pUserCut, isTryUnfold); // solve this cluster and add all unfolded clusters to provided list + // clu.reset(); //empty current cluster + } // digits loop for current chamber + } // chambers loop + return; +} // Dig2Clu() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Clusterer::FormClu(Cluster* pClu, int pDig, std::vector* pDigList, TMatrixF* pDigMap) +{ + // Forms the initial cluster as a combination of all adjascent digits. Starts from the given digit then calls itself recursevly for all neighbours. + // Arguments: pClu - pointer to cluster being formed + // Returns: none + pClu->digAdd(&pDigList->at(pDig)); // take this digit in cluster + int cnt = 0; + int cx[4]; + int cy[4]; + int padChX = 0; + int padChY = 0; + int module = 0; + o2::hmpid::Digit::pad2Absolute(pDigList->at(pDig).getPadID(), &module, &padChX, &padChY); + + if (padChX > Param::kMinPx) { + cx[cnt] = padChX - 1; + cy[cnt] = padChY; + cnt++; + } // left + if (padChX < Param::kMaxPcx) { + cx[cnt] = padChX + 1; + cy[cnt] = padChY; + cnt++; + } // right + if (padChY > Param::kMinPy) { + cx[cnt] = padChX; + cy[cnt] = padChY - 1; + cnt++; + } // down + if (padChY < Param::kMaxPcy) { + cx[cnt] = padChX; + cy[cnt] = padChY + 1; + cnt++; + } // up + + for (int i = 0; i < cnt; i++) { // neighbours loop + pDig = UseDig(cx[i], cy[i], pDigMap); + if (pDig != -1) { + FormClu(pClu, pDig, pDigList, pDigMap); + } // check if this neighbour pad fired and mark it as taken + } // neighbours loop +} // FormClu() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +int Clusterer::UseDig(int padX, int padY, TMatrixF* pPadMap) +{ + // Digit map contains a matrix if digit numbers. + // Main operation in forming initial cluster is done here. Requested digit pointer is returned and this digit marked as taken. + // Arguments: padX,padY - pad number + // pDigLst - list of digits for one sector + // pDigMap - map of those digits + // Returns: index to digit if not yet used or -1 if used + int iDig = (int)(*pPadMap)(padX, padY); + (*pPadMap)(padX, padY) = -1; // take digit number from the map and reset this map cell to -1 + return iDig; +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// bool Clusterer::IsDigSurvive(Int_t *pUserCut, Digit *pDig)const +bool Clusterer::IsDigSurvive(Digit* pDig) const +{ + // Check if the current digit survive to a riapllied sigma cut + // Arguments: pDig pointer to the current digit + // Returns: kTRUE if charge > mean+n*sigma + /*int iCh = pDig->Ch(); + int iDaqSigCut =(int)fDaqSig->At(iCh)->GetUniqueID(); + if(pUserCut[iCh]<=iDaqSigCut) return kTRUE; + TMatrixF *pM = (TMatrixF*)fDaqSig->At(pDig->Ch()); + float sig = (*pM)(pDig->PadChX(),pDig->PadChY()); + //if(pDig->Q()>pUserCut[iCh]*sig) return kTRUE; to be improved + */ + if (pDig->getQ() > 4.) { + return true; + } else { + return false; + } +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +/*void Clusterer::process(std::vector const& digits, std::vector& clusters, MCLabelContainer const* digitMCTruth) { TStopwatch timerProcess; timerProcess.Start(); @@ -44,4 +171,4 @@ void Clusterer::process(std::vector const& digits, std::vector printf("Timing:\n"); printf("Clusterer::process: "); timerProcess.Print(); -} +}*/ diff --git a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx index af4e5cff6e180..1bc1b3354e492 100644 --- a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx +++ b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx @@ -63,6 +63,8 @@ void HMPIDDigitizer::reset() // this will process hits and fill the digit vector with digits which are finalized void HMPIDDigitizer::process(std::vector const& hits, std::vector& digits) { + printf("*******************In digitizer process**********************************"); + for (auto& hit : hits) { int chamber, pc, px, py; float totalQ; diff --git a/Detectors/HMPID/simulation/src/HmpidCoder2.cxx b/Detectors/HMPID/simulation/src/HmpidCoder2.cxx index 4a9d0bbdf8d79..ff4e532ab22b9 100644 --- a/Detectors/HMPID/simulation/src/HmpidCoder2.cxx +++ b/Detectors/HMPID/simulation/src/HmpidCoder2.cxx @@ -247,12 +247,14 @@ void HmpidCoder2::openOutputStream(const std::string& outputFileName, const std: rdh.linkID = ReadOut::LnkId(eq); rdh.endPointID = 0; std::string outfname; - if (fileFor == "link" || fileFor == "cru") { // for the C-RORC detector the link is equivalent of the CRU endpoint + if (fileFor == "link") { outfname = fmt::format("{}_{}_feeid{}.raw", outputFileName, ReadOut::FlpHostName(eq), ReadOut::FeeId(eq)); } else if (fileFor == "flp") { outfname = fmt::format("{}_{}.raw", outputFileName, ReadOut::FlpHostName(eq)); } else if (fileFor == "all") { outfname = fmt::format("{}.raw", outputFileName); + } else if (fileFor == "cru") { + outfname = fmt::format("{}_{}.raw", outputFileName, ReadOut::FlpHostName(eq)); } else { throw std::runtime_error(fmt::format("unknown raw file grouping option {}", fileFor)); } diff --git a/Detectors/HMPID/workflow/CMakeLists.txt b/Detectors/HMPID/workflow/CMakeLists.txt index a65251c1f4539..4599e692330e9 100644 --- a/Detectors/HMPID/workflow/CMakeLists.txt +++ b/Detectors/HMPID/workflow/CMakeLists.txt @@ -13,6 +13,9 @@ o2_add_library(HMPIDWorkflow SOURCES src/DataDecoderSpec.cxx src/DataDecoderSpec2.cxx src/DigitsToRawSpec.cxx + src/DigitsToClustersSpec.cxx + src/DigitsToRootSpec.cxx + src/ClustersToRootSpec.cxx src/DumpDigitsSpec.cxx src/PedestalsCalculationSpec.cxx src/RawToDigitsSpec.cxx @@ -27,17 +30,14 @@ o2_add_library(HMPIDWorkflow O2::DetectorsRaw O2::HMPIDBase O2::DataFormatsHMP + O2::DataFormatsDCS O2::HMPIDSimulation O2::HMPIDReconstruction) - -#o2_add_executable(recoworkflow -# COMPONENT_NAME hmpid -# SOURCES src/HMPIDRecoWorkflow.cxx -# PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + o2_add_executable(entropy-encoder-workflow - COMPONENT_NAME hhmpid - SOURCES src/entropy-encoder-workflow.cxx - PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + COMPONENT_NAME hhmpid + SOURCES src/entropy-encoder-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) o2_add_executable(read-raw-file-stream-workflow COMPONENT_NAME hmpid @@ -58,7 +58,7 @@ o2_add_executable(raw-to-digits-stream-workflow COMPONENT_NAME hmpid SOURCES src/raw-to-digits-stream-workflow.cxx PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) - + o2_add_executable(dump-digits-stream-workflow COMPONENT_NAME hmpid SOURCES src/dump-digits-stream-workflow.cxx @@ -69,6 +69,21 @@ o2_add_executable(digits-to-raw-workflow SOURCES src/digits-to-raw-workflow.cxx PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) +o2_add_executable(digits-to-clusters-workflow + COMPONENT_NAME hmpid + SOURCES src/digits-to-clusters-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + +o2_add_executable(digits-to-root-workflow + COMPONENT_NAME hmpid + SOURCES src/digits-to-root-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + +o2_add_executable(clusters-to-root-workflow + COMPONENT_NAME hmpid + SOURCES src/clusters-to-root-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + o2_add_executable(digits-to-raw-stream-workflow COMPONENT_NAME hmpid SOURCES src/digits-to-raw-stream-workflow.cxx diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h new file mode 100644 index 0000000000000..0e6d8d7ad0e5a --- /dev/null +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h @@ -0,0 +1,63 @@ +// 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. + +/// +/// \file ClustersToRootSpec.h +/// \author Antonio Franco +/// +/// \brief Definition of a data processor to write Root File from Clusters stream +/// + +#ifndef DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_CLUSTERSTOROOTSPEC_H_ +#define DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_CLUSTERSTOROOTSPEC_H_ + +#include "TTree.h" +#include "TFile.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +#include "HMPIDBase/Common.h" +#include "DataFormatsHMP/Cluster.h" +#include "DataFormatsHMP/Trigger.h" +#include "CommonDataFormat/InteractionRecord.h" + +#include "Framework/WorkflowSpec.h" + +namespace o2 +{ +namespace hmpid +{ + +class ClustersToRootTask : public framework::Task +{ + public: + ClustersToRootTask() = default; + ~ClustersToRootTask() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) override; + + private: + ExecutionTimer mExTimer; + std::vector mTriggers; + std::vector mClusters; + TTree* mTheTree; + std::string mOutRootFileName; + TFile* mfileOut; +}; + +o2::framework::DataProcessorSpec getClustersToRootSpec(std::string inputSpec = "HMP/CLUSTERS"); + +} // end namespace hmpid +} // end namespace o2 + +#endif diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h new file mode 100644 index 0000000000000..1f6473c95e7c2 --- /dev/null +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h @@ -0,0 +1,59 @@ +// 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. + +/// +/// \file DatDecoderSpec.h +/// \author Andrea Ferrero +/// +/// \brief Definition of a data processor to run the raw decoding +/// + +#ifndef DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_DIGITSTOCLUSTERSPEC_H_ +#define DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_DIGITSTOCLUSTERSPEC_H_ + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +#include "HMPIDBase/Common.h" +#include "HMPIDReconstruction/Clusterer.h" +#include "Framework/WorkflowSpec.h" + +namespace o2 +{ +namespace hmpid +{ + +class DigitsToClustersTask : public framework::Task +{ + public: + DigitsToClustersTask() = default; + ~DigitsToClustersTask() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) override; + + private: + std::string mSigmaCutPar; + float mSigmaCut[7] = {4, 4, 4, 4, 4, 4, 4}; + + o2::hmpid::Clusterer* mRec; + long mDigitsReceived; + + ExecutionTimer mExTimer; + void strToFloatsSplit(std::string s, std::string delimiter, float* res, int maxElem = 7); +}; + +o2::framework::DataProcessorSpec getDigitsToClustersSpec(std::string inputSpec = "HMP/DIGITS"); + +} // end namespace hmpid +} // end namespace o2 + +#endif diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToRootSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToRootSpec.h new file mode 100644 index 0000000000000..45edadfa409fa --- /dev/null +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToRootSpec.h @@ -0,0 +1,63 @@ +// 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. + +/// +/// \file DigitsToRootSpec.h +/// \author Antonio Franco +/// +/// \brief Definition of a data processor to store digits in Root file +/// + +#ifndef DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_DIGITSTOROOTSPEC_H_ +#define DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_DIGITSTOROOTSPEC_H_ + +#include "TTree.h" +#include "TFile.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +#include "HMPIDBase/Common.h" +#include "DataFormatsHMP/Digit.h" +#include "DataFormatsHMP/Trigger.h" +#include "CommonDataFormat/InteractionRecord.h" + +#include "Framework/WorkflowSpec.h" + +namespace o2 +{ +namespace hmpid +{ + +class DigitsToRootTask : public framework::Task +{ + public: + DigitsToRootTask() = default; + ~DigitsToRootTask() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) override; + + private: + ExecutionTimer mExTimer; + std::vector mTriggers; + std::vector mDigits; + TFile* mfileOut; + TTree* mTheTree; + std::string mOutRootFileName; +}; + +o2::framework::DataProcessorSpec getDigitsToRootSpec(std::string inputSpec = "HMP/DIGITS"); + +} // end namespace hmpid +} // end namespace o2 + +#endif diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/PedestalsCalculationSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/PedestalsCalculationSpec.h index dc0c952484dd6..296de25a22e25 100644 --- a/Detectors/HMPID/workflow/include/HMPIDWorkflow/PedestalsCalculationSpec.h +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/PedestalsCalculationSpec.h @@ -33,31 +33,31 @@ class PedestalsCalculationTask : public framework::Task void init(framework::InitContext& ic) final; void run(framework::ProcessingContext& pc) final; void decodeTF(framework::ProcessingContext& pc); - void decodeReadout(framework::ProcessingContext& pc); - void decodeRawFile(framework::ProcessingContext& pc); void endOfStream(framework::EndOfStreamContext& ec) override; private: void recordPedInCcdb(); + void recordPedInFiles(); + void recordPedInDcsCcdb(); private: HmpidDecoder2* mDeco; long mTotalDigits; long mTotalFrames; - std::string mPedestalsBasePath; float mSigmaCut; std::string mPedestalTag; - + bool mWriteToFiles; + std::string mPedestalsBasePath; o2::ccdb::CcdbApi mDBapi; std::map mDbMetadata; // can be empty + std::string mPedestalsCCDBBasePath; bool mWriteToDB; bool mFastAlgorithm; - ExecutionTimer mExTimer; }; o2::framework::DataProcessorSpec getPedestalsCalculationSpec(std::string inputSpec = "TF:HMP/RAWDATA"); -//o2::framework::DataProcessorSpec getDecodingSpec(); + } // end namespace hmpid } // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx b/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx new file mode 100644 index 0000000000000..6b39d74b41f87 --- /dev/null +++ b/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx @@ -0,0 +1,153 @@ +// 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. + +/// \file ClustersToRootSpec.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 20 nov 2021 +/// \brief Implementation of a data processor to produce Root File from Clusters Stream +/// + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DPLUtils/DPLRawParser.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" + +#include "TTree.h" +#include "TFile.h" + +#include + +#include "Framework/DataRefUtils.h" +#include "Framework/InputSpec.h" +#include "Framework/CallbackService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" +#include "Framework/Logger.h" +#include "Framework/InputRecordWalker.h" + +#include "Headers/RAWDataHeader.h" +#include "DetectorsRaw/RDHUtils.h" + +#include "HMPIDBase/Geo.h" +#include "HMPIDWorkflow/ClustersToRootSpec.h" + +namespace o2 +{ +namespace hmpid +{ + +using namespace o2; +using namespace o2::framework; +using RDH = o2::header::RDHAny; + +//======================= +// Data decoder +void ClustersToRootTask::init(framework::InitContext& ic) +{ + LOG(info) << "[HMPID Write Root File From Clusters stream - init()]"; + + // get line command options + mOutRootFileName = ic.options().get("out-file"); + + mTriggers.clear(); + mClusters.clear(); + + TString filename; + TString tit; + + filename = TString::Format("%s", mOutRootFileName.c_str()); + LOG(info) << "Create the ROOT file " << filename.Data(); + mfileOut = new TFile(TString::Format("%s", filename.Data()), "RECREATE"); + tit = TString::Format("HMPID Clusters File Decoding"); + mTheTree = new TTree("o2hmp", tit); + mTheTree->Branch("InteractionRecords", &mTriggers); + mTheTree->Branch("HMPIDClusters", &mClusters); + + mExTimer.start(); + return; +} + +void ClustersToRootTask::run(framework::ProcessingContext& pc) +{ + std::vector triggers; + std::vector clusters; + + for (auto const& ref : InputRecordWalker(pc.inputs())) { + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "INTRECORDS1"}})) { + triggers = pc.inputs().get>(ref); + // LOG(info) << "We receive triggers =" << triggers.size(); + } + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "CLUSTERS"}})) { + clusters = pc.inputs().get>(ref); + // LOG(info) << "The size of the vector =" << clusters.size(); + } + + for (int i = 0; i < triggers.size(); i++) { + // LOG(info) << "Trigger Event Orbit = " << triggers[i].getOrbit() << " BC = " << triggers[i].getBc(); + int startClusterIndex = mClusters.size(); + int numberOfClusters = 0; + for (int j = triggers[i].getFirstEntry(); j <= triggers[i].getLastEntry(); j++) { + mClusters.push_back(clusters[j]); // append the cluster + numberOfClusters++; + } + mTriggers.push_back(triggers[i]); + mTriggers.back().setDataRange(startClusterIndex, numberOfClusters); + } + } + mExTimer.stop(); + return; +} + +void ClustersToRootTask::endOfStream(framework::EndOfStreamContext& ec) +{ + mExTimer.logMes("Received an End Of Stream !"); + for (int i = 0; i < mClusters.size(); i++) { + mClusters.at(i).print(); + } + mTheTree->Fill(); + mTheTree->Write(); + mfileOut->Close(); + mExTimer.logMes("Register Tree ! "); + return; +} + +//_________________________________________________________________________________________________ +o2::framework::DataProcessorSpec getClustersToRootSpec(std::string inputSpec) +{ + std::vector inputs; + inputs.emplace_back("clusters", o2::header::gDataOriginHMP, "CLUSTERS", 0, Lifetime::Timeframe); + inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS1", 0, Lifetime::Timeframe); + + std::vector outputs; + + return DataProcessorSpec{ + "HMP-ClustersToRoot", + inputs, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{{"out-file", VariantType::String, "hmpClusters.root", {"name of the output file"}}}}; +} + +} // namespace hmpid +} // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/DataDecoderSpec2.cxx b/Detectors/HMPID/workflow/src/DataDecoderSpec2.cxx index 123c5955d328b..6c71bb9eaf482 100644 --- a/Detectors/HMPID/workflow/src/DataDecoderSpec2.cxx +++ b/Detectors/HMPID/workflow/src/DataDecoderSpec2.cxx @@ -197,7 +197,7 @@ void DataDecoderTask2::decodeTF(framework::ProcessingContext& pc) } DPLRawParser parser(inputs, o2::framework::select("TF:HMP/RAWDATA")); - //mDeco->mDigits.clear(); + // mDeco->mDigits.clear(); for (auto it = parser.begin(), end = parser.end(); it != end; ++it) { int pointerToTheFirst = mDeco->mDigits.size(); uint32_t* theBuffer = (uint32_t*)it.raw(); diff --git a/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx b/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx new file mode 100644 index 0000000000000..0481b89e04797 --- /dev/null +++ b/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx @@ -0,0 +1,171 @@ +// 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. + +/// \file DigitsToRootSpec.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 01 feb 2021 +/// \brief Implementation of a data processor to Clusterize the Digits +/// + +#include +#include +#include +#include +#include +#include +#include + +#include "HMPIDWorkflow/DigitsToClustersSpec.h" +#include "Framework/CallbackService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" +#include "Framework/Logger.h" +#include "Framework/DataRefUtils.h" +#include "Framework/InputRecordWalker.h" + +#include "Headers/RAWDataHeader.h" +#include "DetectorsRaw/RDHUtils.h" +#include "DPLUtils/DPLRawParser.h" + +#include "DataFormatsHMP/Digit.h" +#include "DataFormatsHMP/Trigger.h" +#include "DataFormatsHMP/Cluster.h" +#include "HMPIDBase/Geo.h" + +namespace o2 +{ +namespace hmpid +{ + +using namespace o2; +using namespace o2::header; +using namespace o2::framework; +using RDH = o2::header::RDHAny; + +// Splits a string in float array for string delimiter, TODO: Move this in a HMPID common library +void DigitsToClustersTask::strToFloatsSplit(std::string s, std::string delimiter, float* res, int maxElem) +{ + int index = 0; + size_t pos_start = 0; + size_t pos_end; + size_t delim_len = delimiter.length(); + std::string token; + while ((pos_end = s.find(delimiter, pos_start)) != std::string::npos) { + token = s.substr(pos_start, pos_end - pos_start); + pos_start = pos_end + delim_len; + res[index++] = std::stof(token); + if (index == maxElem) { + return; + } + } + res[index++] = (std::stof(s.substr(pos_start))); + return; +} + +//======================= +// +void DigitsToClustersTask::init(framework::InitContext& ic) +{ + LOG(info) << "[HMPID Clusterization - init() ] "; + mSigmaCutPar = ic.options().get("sigma-cut"); + if (mSigmaCutPar != "") { + strToFloatsSplit(mSigmaCutPar, ",", mSigmaCut, 7); + } + + mDigitsReceived = 0; + mRec = new o2::hmpid::Clusterer(); + mExTimer.start(); + return; +} + +void DigitsToClustersTask::run(framework::ProcessingContext& pc) +{ + LOG(info) << "[HMPID DClusterization - run() ] Enter ..."; + std::vector triggers; + std::vector digits; + std::vector oneEventDigits; + std::vector oneEventClusters; + std::vector clusters; + std::vector clustersTrigger; + + for (auto const& ref : InputRecordWalker(pc.inputs())) { + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{gDataOriginHMP, "INTRECORDS"}})) { + triggers = pc.inputs().get>(ref); + LOG(info) << "We receive triggers =" << triggers.size(); + } + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{gDataOriginHMP, "DIGITS"}})) { + digits = pc.inputs().get>(ref); + LOG(info) << "The size of the vector =" << digits.size(); + mDigitsReceived += digits.size(); + } + + clustersTrigger.clear(); + for (int i = 0; i < triggers.size(); i++) { + oneEventDigits.clear(); + LOG(info) << "Trigger Event Orbit = " << triggers[i].getOrbit() << " BC = " << triggers[i].getBc(); + for (int j = triggers[i].getFirstEntry(); j <= triggers[i].getLastEntry(); j++) { + oneEventDigits.push_back(digits[j]); + } + // ------ Now we have the Digit Array for one Trigger event + oneEventClusters.clear(); + mRec->Dig2Clu(&oneEventDigits, &oneEventClusters, mSigmaCut, true); + if (oneEventClusters.size() > 0) { // we have something to store + clustersTrigger.push_back(triggers[i]); // copy the Interaction record in the new triggers vector + clustersTrigger.back().setDataRange(clusters.size(), oneEventClusters.size()); // set the data range to the chunch of clusters + for (int j = 0; j < oneEventClusters.size(); j++) { + clusters.push_back(oneEventClusters[j]); + } // append clusters + } + LOG(info) << "Clusters for event = " << oneEventClusters.size(); + } + + pc.outputs().snapshot(o2::framework::Output{"HMP", "CLUSTERS", 0, o2::framework::Lifetime::Timeframe}, clusters); + pc.outputs().snapshot(o2::framework::Output{"HMP", "INTRECORDS1", 0, o2::framework::Lifetime::Timeframe}, clustersTrigger); + } + + mExTimer.elapseMes("Clusterization of Digits received = " + std::to_string(mDigitsReceived)); + return; +} + +void DigitsToClustersTask::endOfStream(framework::EndOfStreamContext& ec) +{ + mExTimer.stop(); + mExTimer.logMes("End Clusterization ! digits = " + std::to_string(mDigitsReceived)); + return; +} + +//_________________________________________________________________________________________________ +o2::framework::DataProcessorSpec getDigitsToClustersSpec(std::string inputSpec) +{ + std::vector inputs; + inputs.emplace_back("digits", o2::header::gDataOriginHMP, "DIGITS", 0, Lifetime::Timeframe); + inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS", 0, Lifetime::Timeframe); + + std::vector outputs; + outputs.emplace_back("HMP", "CLUSTERS", 0, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("HMP", "INTRECORDS1", 0, o2::framework::Lifetime::Timeframe); + + return DataProcessorSpec{ + "HMP-Clusterization", + inputs, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{{"sigma-cut", VariantType::String, "", {"sigmas as comma separated list"}}}}; +} + +} // namespace hmpid +} // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/DigitsToRootSpec.cxx b/Detectors/HMPID/workflow/src/DigitsToRootSpec.cxx new file mode 100644 index 0000000000000..f9a6b7d7038b0 --- /dev/null +++ b/Detectors/HMPID/workflow/src/DigitsToRootSpec.cxx @@ -0,0 +1,149 @@ +// 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. + +/// \file DigitsToRootSpec.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 20 nov 2021 +/// \brief Implementation of a data processor to produce Root File from Digits stream +/// + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DPLUtils/DPLRawParser.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" + +#include "TTree.h" +#include "TFile.h" + +#include + +#include "Framework/DataRefUtils.h" +#include "Framework/InputSpec.h" +#include "Framework/CallbackService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" +#include "Framework/Logger.h" +#include "Framework/InputRecordWalker.h" + +#include "Headers/RAWDataHeader.h" +#include "DetectorsRaw/RDHUtils.h" + +#include "HMPIDBase/Geo.h" +#include "HMPIDWorkflow/DigitsToRootSpec.h" + +namespace o2 +{ +namespace hmpid +{ + +using namespace o2; +using namespace o2::framework; +using RDH = o2::header::RDHAny; + +//======================= +void DigitsToRootTask::init(framework::InitContext& ic) +{ + LOG(info) << "[HMPID Write Root File From Digits stream - init()]"; + + // get line command options + mOutRootFileName = ic.options().get("out-file"); + + mTriggers.clear(); + mDigits.clear(); + + TString filename = TString::Format("%s", mOutRootFileName.c_str()); + TString tit = TString::Format("HMPID Digits File Decoding"); + + LOG(info) << "Create the ROOT file " << filename.Data(); + mfileOut = new TFile(TString::Format("%s", filename.Data()), "RECREATE"); + + mTheTree = new TTree("o2hmp", tit); + mTheTree->Branch("InteractionRecords", &mTriggers); + mTheTree->Branch("HMPIDDigits", &mDigits); + + mExTimer.start(); + return; +} + +void DigitsToRootTask::run(framework::ProcessingContext& pc) +{ + std::vector triggers; + std::vector digits; + + for (auto const& ref : InputRecordWalker(pc.inputs())) { + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "INTRECORDS"}})) { + triggers = pc.inputs().get>(ref); + LOG(info) << "We receive triggers =" << triggers.size(); + } + if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "DIGITS"}})) { + digits = pc.inputs().get>(ref); + LOG(info) << "The size of the vector =" << digits.size(); + } + + for (int i = 0; i < triggers.size(); i++) { + LOG(info) << "Trigger Event Orbit = " << triggers[i].getOrbit() << " BC = " << triggers[i].getBc(); + int startDigitsIndex = mDigits.size(); + int numberOfDigits = 0; + for (int j = triggers[i].getFirstEntry(); j <= triggers[i].getLastEntry(); j++) { + mDigits.push_back(digits[j]); // append the cluster + numberOfDigits++; + } + mTriggers.push_back(triggers[i]); + mTriggers.back().setDataRange(startDigitsIndex, numberOfDigits); + } + } + mExTimer.stop(); + return; +} + +void DigitsToRootTask::endOfStream(framework::EndOfStreamContext& ec) +{ + mExTimer.logMes("Received an End Of Stream !"); + LOG(info) << "The size of digits vector =" << mDigits.size(); + mTheTree->Fill(); + mTheTree->Write(); + mfileOut->Close(); + mExTimer.logMes("Register Tree ! "); + return; +} + +//_________________________________________________________________________________________________ +o2::framework::DataProcessorSpec getDigitsToRootSpec(std::string inputSpec) +{ + std::vector inputs; + inputs.emplace_back("clusters", o2::header::gDataOriginHMP, "DIGITS", 0, Lifetime::Timeframe); + inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS", 0, Lifetime::Timeframe); + + std::vector outputs; + + return DataProcessorSpec{ + "HMP-DigitsToRoot", + inputs, + outputs, + AlgorithmSpec{adaptFromTask()}, + Options{{"out-file", VariantType::String, "hmpDigits.root", {"name of the output file"}}}}; +} + +} // namespace hmpid +} // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/PedestalsCalculationSpec.cxx b/Detectors/HMPID/workflow/src/PedestalsCalculationSpec.cxx index fac947ff69f95..3f1464af63b0f 100644 --- a/Detectors/HMPID/workflow/src/PedestalsCalculationSpec.cxx +++ b/Detectors/HMPID/workflow/src/PedestalsCalculationSpec.cxx @@ -25,6 +25,7 @@ #include #include "TTree.h" +#include "TString.h" #include "TFile.h" #include "TMath.h" #include "TMatrixF.h" @@ -46,14 +47,15 @@ #include "Headers/RAWDataHeader.h" #include "DetectorsRaw/RDHUtils.h" #include "DPLUtils/DPLRawParser.h" +#include "CommonUtils/NameConf.h" #include "DataFormatsHMP/Digit.h" #include "DataFormatsHMP/Trigger.h" #include "HMPIDBase/Geo.h" #include "HMPIDReconstruction/HmpidDecoder2.h" #include "HMPIDWorkflow/PedestalsCalculationSpec.h" -#include "CommonUtils/NameConf.h" +#include "DataFormatsDCS/DCSConfigObject.h" namespace o2 { namespace hmpid @@ -64,7 +66,7 @@ using namespace o2::framework; using RDH = o2::header::RDHAny; //======================= -// Data decoder +// Init void PedestalsCalculationTask::init(framework::InitContext& ic) { @@ -76,28 +78,43 @@ void PedestalsCalculationTask::init(framework::InitContext& ic) mTotalFrames = 0; mSigmaCut = ic.options().get("sigmacut"); + mWriteToFiles = ic.options().get("use-files"); mPedestalsBasePath = ic.options().get("files-basepath"); + mPedestalsCCDBBasePath = mPedestalsBasePath; + + mWriteToDB = ic.options().get("use-ccdb"); + if (mWriteToDB) { + mDBapi.init(ic.options().get("ccdb-uri")); // or http://localhost:8080 for a local installation + mWriteToDB = mDBapi.isHostReachable() ? true : false; + } mPedestalTag = ic.options().get("pedestals-tag"); - mDBapi.init(ic.options().get("ccdb-uri")); // or http://localhost:8080 for a local installation - mWriteToDB = mDBapi.isHostReachable() ? true : false; mFastAlgorithm = ic.options().get("fast-decode"); mExTimer.start(); + LOG(info) << "Calculate Ped/Thresh." + (mWriteToDB ? " Store in DCSCCDB at " + mPedestalsBasePath + " with Tag:" + mPedestalTag : " CCDB not used !"); return; } void PedestalsCalculationTask::run(framework::ProcessingContext& pc) { decodeTF(pc); - // TODO: accept other types of Raw Streams ... - // decodeReadout(pc); - // decodeRawFile(pc);ccdb - mExTimer.elapseMes("Decoding... Digits decoded = " + std::to_string(mTotalDigits) + " Frames received = " + std::to_string(mTotalFrames)); return; } void PedestalsCalculationTask::endOfStream(framework::EndOfStreamContext& ec) +{ + if (mWriteToDB) { + recordPedInDcsCcdb(); + } + if (mWriteToFiles) { + recordPedInFiles(); + } + mExTimer.stop(); + return; +} + +void PedestalsCalculationTask::recordPedInFiles() { double Average; double Variance; @@ -116,7 +133,11 @@ void PedestalsCalculationTask::endOfStream(framework::EndOfStreamContext& ec) } sprintf(padsFileName, "%s_%d.dat", mPedestalsBasePath.c_str(), e); FILE* fpads = fopen(padsFileName, "w"); - // TODO: Add controls on the file open + if (fpads == nullptr) { + mExTimer.logMes("error creating the file = " + std::string(padsFileName)); + LOG(error) << "error creating the file = " << padsFileName; + return; + } for (int c = 0; c < Geo::N_COLUMNS; c++) { for (int d = 0; d < Geo::N_DILOGICS; d++) { for (int h = 0; h < Geo::N_CHANNELS; h++) { @@ -146,18 +167,92 @@ void PedestalsCalculationTask::endOfStream(framework::EndOfStreamContext& ec) fclose(fpads); } mExTimer.logMes("End Writing the pedestals ! Digits decoded = " + std::to_string(mTotalDigits) + " Frames received = " + std::to_string(mTotalFrames)); + return; +} - if (mWriteToDB) { - recordPedInCcdb(); +void PedestalsCalculationTask::recordPedInDcsCcdb() +{ + // create the root structure + LOG(info) << "Store Pedestals in DCS CCDB "; + + float xb, yb, ch, Samples; + double SumOfCharge, SumOfSquares, Average, Variance; + uint32_t Pedestal, Threshold, PedThr; + + o2::dcs::DCSconfigObject_t pedestalsConfig; + +// Setup dimensions for Equipment granularity with 48 channels/dilogic +#define PEDTHFORMAT "%05X," +#define COLUMNTAIL false +#define PEDTHBYTES 6 + int bufferDim = PEDTHBYTES * Geo::N_CHANNELS * Geo::N_DILOGICS * Geo::N_COLUMNS + 10; + char* outBuffer = new char[bufferDim]; + char* inserPtr; + + for (int e = 0; e < Geo::MAXEQUIPMENTS; e++) { + if (mDeco->getAverageEventSize(e) == 0) { // skip the empty equipment + continue; + } + inserPtr = outBuffer; + // algoritm based on equipment granularity + for (int c = 0; c < Geo::N_COLUMNS; c++) { + for (int d = 0; d < Geo::N_DILOGICS; d++) { + for (int h = 0; h < Geo::N_CHANNELS; h++) { + Samples = (double)mDeco->getChannelSamples(e, c, d, h); + SumOfCharge = mDeco->getChannelSum(e, c, d, h); + SumOfSquares = mDeco->getChannelSquare(e, c, d, h); + if (Samples > 0) { + Average = SumOfCharge / Samples; + Variance = sqrt(abs((Samples * SumOfSquares) - (SumOfCharge * SumOfCharge))) / Samples; + } else { + Average = 0; + Variance = 0; + } + Pedestal = (uint32_t)Average; + Threshold = (uint32_t)(Variance * mSigmaCut + Average); + PedThr = ((Threshold & 0x001FF) << 9) | (Pedestal & 0x001FF); + sprintf(inserPtr, PEDTHFORMAT, PedThr); + inserPtr += PEDTHBYTES; + } + if (COLUMNTAIL) { + for (int h = 48; h < 64; h++) { + sprintf(inserPtr, PEDTHFORMAT, 0); + inserPtr += PEDTHBYTES; + } + } + } + } + mExTimer.logMes("End write the equipment = " + std::to_string(e)); + sprintf(inserPtr, "%05X\n", 0xA0A0A); // The closure value + inserPtr += 6; + *inserPtr = '\0'; // close the string rap. + o2::dcs::addConfigItem(pedestalsConfig, "Equipment" + std::to_string(e), (const char*)outBuffer); } - mExTimer.stop(); + + struct timeval tp; + gettimeofday(&tp, nullptr); + uint64_t ms = tp.tv_sec * 1000 + tp.tv_usec / 1000; + uint64_t minTimeStamp = ms; + uint64_t maxTimeStamp = ms + 1; + + char filename[1024]; + sprintf(filename, "%s/%s/PedThre.root", mPedestalsCCDBBasePath.c_str(), mPedestalTag.c_str()); + mExTimer.logMes("File name = >" + std::string(filename) + "< (" + mPedestalsCCDBBasePath + "," + mPedestalTag); + TFile outputFile(filename, "recreate"); + outputFile.WriteObjectAny(&pedestalsConfig, "std::vector", "DCSConfig"); + outputFile.Close(); + + mDbMetadata.emplace("Tag", mPedestalTag.c_str()); + mDBapi.storeAsTFileAny(&pedestalsConfig, filename, mDbMetadata, minTimeStamp, maxTimeStamp); + + mExTimer.logMes("End Writing the pedestals ! Digits decoded = " + std::to_string(mTotalDigits) + " Frames received = " + std::to_string(mTotalFrames)); + return; } void PedestalsCalculationTask::recordPedInCcdb() { // create the root structure - LOG(info) << "Store Pedestals in ccdb "; float xb, yb, ch; @@ -167,8 +262,8 @@ void PedestalsCalculationTask::recordPedInCcdb() TObjArray aPedestals(Geo::N_MODULES); for (int i = 0; i < Geo::N_MODULES; i++) { - aSigmas.AddAt(new TMatrixF(160, 144), i); - aPedestals.AddAt(new TMatrixF(160, 144), i); + aSigmas.AddAt(new TMatrixF(Geo::N_XROWS, Geo::N_YCOLS), i); + aPedestals.AddAt(new TMatrixF(Geo::N_XROWS, Geo::N_YCOLS), i); } for (int m = 0; m < o2::hmpid::Geo::N_MODULES; m++) { @@ -197,7 +292,6 @@ void PedestalsCalculationTask::recordPedInCcdb() struct timeval tp; gettimeofday(&tp, nullptr); uint64_t ms = tp.tv_sec * 1000 + tp.tv_usec / 1000; - uint64_t minTimeStamp = ms; uint64_t maxTimeStamp = ms + 1; @@ -205,7 +299,7 @@ void PedestalsCalculationTask::recordPedInCcdb() if (mDeco->getAverageEventSize(i * 2) == 0 && mDeco->getAverageEventSize(i * 2 + 1) == 0) { continue; // If no events skip the chamber } - TString filename = TString::Format("HMP/Pedestals/%s/Mean_%d", mPedestalTag.c_str(), i); + TString filename = TString::Format("%s/%s/Mean_%d", mPedestalsCCDBBasePath.c_str(), mPedestalTag.c_str(), i); mDbMetadata.emplace("Tag", mPedestalTag.c_str()); mDBapi.storeAsTFileAny(aPedestals.At(i), filename.Data(), mDbMetadata, minTimeStamp, maxTimeStamp); } @@ -213,7 +307,7 @@ void PedestalsCalculationTask::recordPedInCcdb() if (mDeco->getAverageEventSize(i * 2) == 0 && mDeco->getAverageEventSize(i * 2 + 1) == 0) { continue; // If no events skip the chamber } - TString filename = TString::Format("HMP/Pedestals/%s/Sigma_%d", mPedestalTag.c_str(), i); + TString filename = TString::Format("%s/%s/Sigma_%d", mPedestalsCCDBBasePath.c_str(), mPedestalTag.c_str(), i); mDbMetadata.emplace("Tag", mPedestalTag.c_str()); mDBapi.storeAsTFileAny(aSigmas.At(i), filename.Data(), mDbMetadata, minTimeStamp, maxTimeStamp); } @@ -225,7 +319,6 @@ void PedestalsCalculationTask::recordPedInCcdb() void PedestalsCalculationTask::decodeTF(framework::ProcessingContext& pc) { LOG(debug) << "*********** In decodeTF **************"; - // get the input buffer auto& inputs = pc.inputs(); DPLRawParser parser(inputs, o2::framework::select("TF:HMP/RAWDATA")); @@ -240,79 +333,13 @@ void PedestalsCalculationTask::decodeTF(framework::ProcessingContext& pc) } } catch (int e) { // The stream end ! - LOG(debug) << "End Fast Page decoding !"; + LOG(debug) << "End Page decoding !"; } mTotalFrames++; mTotalDigits += mDeco->mDigits.size(); } return; } -//pc.outputs().make -//_________________________________________________________________________________________________ -// the decodeReadout() function processes the messages generated by o2-mch-cru-page-reader-workflow -void PedestalsCalculationTask::decodeReadout(framework::ProcessingContext& pc) -{ - LOG(info) << "*********** In decode readout **************"; - - // get the input buffer - auto& inputs = pc.inputs(); - DPLRawParser parser(inputs, o2::framework::select("readout:HMP/RAWDATA")); - // DPLRawParser parser(inputs, o2::framework::select("HMP/readout")); - - for (auto it = parser.begin(), end = parser.end(); it != end; ++it) { - uint32_t* theBuffer = (uint32_t*)it.raw(); - mDeco->setUpStream(theBuffer, it.size() + it.offset()); - try { - if (mFastAlgorithm) { - mDeco->decodePageFast(&theBuffer); - } else { - mDeco->decodePage(&theBuffer); - } - } catch (int e) { - // The stream end ! - LOG(debug) << "End Fast Page decoding !"; - } - } - return; -} - -// the decodeReadout() function processes the messages generated by o2-mch-cru-page-reader-workflow -void PedestalsCalculationTask::decodeRawFile(framework::ProcessingContext& pc) -{ - LOG(info) << "*********** In decode rawfile **************"; - - for (auto&& input : pc.inputs()) { - if (input.spec->binding == "file") { - const header::DataHeader* header = o2::header::get(input.header); - if (!header) { - return; - } - - auto const* raw = input.payload; - size_t payloadSize = header->payloadSize; - - LOG(info) << " payloadSize=" << payloadSize; - if (payloadSize == 0) { - return; - } - - uint32_t* theBuffer = (uint32_t*)input.payload; - int pagesize = header->payloadSize; - mDeco->setUpStream(theBuffer, pagesize); - try { - if (mFastAlgorithm) { - mDeco->decodePageFast(&theBuffer); - } else { - mDeco->decodePage(&theBuffer); - } - } catch (int e) { - // The stream end ! - LOG(debug) << "End Fast Page decoding !"; - } - } - } - return; -} //_________________________________________________________________________________________________ o2::framework::DataProcessorSpec getPedestalsCalculationSpec(std::string inputSpec) @@ -320,23 +347,19 @@ o2::framework::DataProcessorSpec getPedestalsCalculationSpec(std::string inputSp std::vector inputs; inputs.emplace_back("TF", o2::framework::ConcreteDataTypeMatcher{"HMP", "RAWDATA"}, o2::framework::Lifetime::Timeframe); - inputs.emplace_back("file", o2::framework::ConcreteDataTypeMatcher{"ROUT", "RAWDATA"}, o2::framework::Lifetime::Timeframe); - // inputs.emplace_back("readout", o2::header::gDataOriginHMP, "RAWDATA", 0, Lifetime::Timeframe); - // inputs.emplace_back("readout", o2::header::gDataOriginHMP, "RAWDATA", 0, Lifetime::Timeframe); - // inputs.emplace_back("rawfile", o2::header::gDataOriginHMP, "RAWDATA", 0, Lifetime::Timeframe); std::vector outputs; - // outputs.emplace_back("HMP", "DIGITS", 0, o2::framework::Lifetime::Timeframe); return DataProcessorSpec{ - "HMP-DataDecoder", + "HMP-PestalsCalculation", o2::framework::select(inputSpec.c_str()), outputs, AlgorithmSpec{adaptFromTask()}, - Options{{"files-basepath", VariantType::String, "/tmp/hmpPedThr", {"Name of the Base Path of Pedestals/Thresholds files."}}, - {"use-ccdb", VariantType::Bool, false, {"Register the Pedestals/Threshold values into the CCDB"}}, - {"ccdb-uri", VariantType::String, o2::base::NameConf::getCCDBServer(), {"URI for the CCDB access."}}, - {"fast-decode", VariantType::Bool, false, {"Use the fast algorithm. (error 0.8%)"}}, + Options{{"files-basepath", VariantType::String, "HMP/Config", {"Name of the Base Path of Pedestals/Thresholds files."}}, + {"use-files", VariantType::Bool, false, {"Register the Pedestals/Threshold values into ASCII files"}}, + {"use-ccdb", VariantType::Bool, true, {"Register the Pedestals/Threshold values into the CCDB"}}, + {"ccdb-uri", VariantType::String, "http://ccdb-test.cern.ch:8080", {"URI for the CCDB access."}}, + {"fast-decode", VariantType::Bool, true, {"Use the fast algorithm. (error 0.8%)"}}, {"pedestals-tag", VariantType::String, "Latest", {"The tag applied to this set of pedestals/threshold values"}}, {"sigmacut", VariantType::Float, 4.0f, {"Sigma values for the Thresholds calculation."}}}}; } diff --git a/Detectors/HMPID/workflow/src/RawToDigitsSpec.cxx b/Detectors/HMPID/workflow/src/RawToDigitsSpec.cxx index 47f95b04e954b..d94085917db73 100644 --- a/Detectors/HMPID/workflow/src/RawToDigitsSpec.cxx +++ b/Detectors/HMPID/workflow/src/RawToDigitsSpec.cxx @@ -235,51 +235,6 @@ void RawToDigitsTask::writeResults() sort(mEvents.begin(), mEvents.end()); // sort the events mExTimer.logMes("Sorted Events = " + std::to_string(mEvents.size())); - /* ------ ROOT file version 1 ---------- - o2::hmpid::Digit digit; - o2::hmpid::Trigger event; - TString filename; - TString tit; - - filename = TString::Format("%s", mOutRootFileName.c_str()); - LOG(info) << "Create the ROOT file " << filename.Data(); - TFile mfileOut(TString::Format("%s", filename.Data()), "RECREATE"); - TTree* theTree; - TBranch* theDigits; - TBranch* theEvents; - tit = TString::Format("HMPID Raw File Decoding"); - theTree = new TTree("o2hmp", tit); - - theDigits = theTree->Branch("HMPDigit", &digit, sizeof(o2::hmpid::Digit), 1); - theEvents = theTree->Branch("InteractionRecords", &event, sizeof(o2::hmpid::Trigger), 1); - - o2::hmpid::Trigger prevEvent = mEvents[0]; - uint32_t theFirstDigit = 0; - uint32_t theLastDigit = 0; - for (int e = 0; e < mEvents.size(); e++) { - LOG(info) << "Manage event " << mEvents[e]; - if (prevEvent != mEvents[e]) { // changes the event Flush It - event = prevEvent; - event.setDataRange(theFirstDigit, theLastDigit-theFirstDigit); - theEvents->Fill(); - theFirstDigit = theLastDigit; - prevEvent = mEvents[e]; - } - int first = mEvents[e].getFirstEntry(); - int last = mEvents[e].getLastEntry(); - for(int idx = first; idx <= last; idx++) { - digit = mAccumulateDigits[idx]; - theDigits->Fill(); - theLastDigit++; - } - } - event = prevEvent; - event.setDataRange(theFirstDigit, theLastDigit-theFirstDigit); - theEvents->Fill(); - theTree->Write(); - mfileOut.Close(); - -------------------------- */ - // ---------- ROOT file version 2 --------------- TString filename; TString tit; diff --git a/Detectors/HMPID/workflow/src/ReadRawFileSpec.cxx b/Detectors/HMPID/workflow/src/ReadRawFileSpec.cxx index 6931a5ac0cd5d..5af790354383a 100644 --- a/Detectors/HMPID/workflow/src/ReadRawFileSpec.cxx +++ b/Detectors/HMPID/workflow/src/ReadRawFileSpec.cxx @@ -79,7 +79,7 @@ void RawFileReaderTask::run(framework::ProcessingContext& pc) sleep(1); while (true) { - //usleep(100); + // usleep(100); mInputFile.read((char*)(&rdh), sizeof(RDH)); // read the next RDH, stop if no more data is available if (mInputFile.fail()) { free(outBuffer); @@ -121,7 +121,7 @@ void RawFileReaderTask::run(framework::ProcessingContext& pc) } bufSize = frameSize; // Set the buffer pointer pc.outputs().snapshot(Output{"HMP", "RAWDATA"}, outBuffer, bufSize); - //std::cout << mExTimer.mTimer.CpuTime() << " " << mExTimer.mLastLogTime << std::endl; + // std::cout << mExTimer.mTimer.CpuTime() << " " << mExTimer.mLastLogTime << std::endl; mExTimer.elapseMes("... Reading... Number of Pages = " + std::to_string(numberOfFrames)); } // while (true) diff --git a/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx b/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx new file mode 100644 index 0000000000000..d150e0a3b1778 --- /dev/null +++ b/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx @@ -0,0 +1,71 @@ +// 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. + +/// \file clusters-to-root-workflow.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 22 nov 2021 +/// + +#include "Framework/WorkflowSpec.h" +#include "Framework/DataSpecUtils.h" +#include "Framework/CallbackService.h" +#include "Framework/ControlService.h" +#include "Framework/Task.h" +#include "Framework/CompletionPolicy.h" +#include "Framework/CompletionPolicyHelpers.h" +#include "Framework/DispatchPolicy.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/Variant.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/NameConf.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "Framework/CallbacksPolicy.h" + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// customize the completion policy +void customize(std::vector& policies) +{ + using o2::framework::CompletionPolicy; + using o2::framework::CompletionPolicyHelpers; + policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName("clusters-hmpid-root", CompletionPolicy::CompletionOp::Consume)); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + std::string keyvaluehelp("Semicolon separated key=value strings ..."); + workflowOptions.push_back(o2::framework::ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {keyvaluehelp}}); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +#include "Framework/runDataProcessing.h" +#include "HMPIDWorkflow/ClustersToRootSpec.h" + +using namespace o2; +using namespace o2::framework; + +WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) +{ + WorkflowSpec specs; + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + DataProcessorSpec consumer = o2::hmpid::getClustersToRootSpec(); + specs.push_back(consumer); + + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + + return specs; +} diff --git a/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx b/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx new file mode 100644 index 0000000000000..513b6adcae70f --- /dev/null +++ b/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx @@ -0,0 +1,67 @@ +// 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. + +/// \file digits-to-cluster-workflow.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 22 nov 2021 +/// + +#include "Framework/WorkflowSpec.h" +#include "Framework/DataSpecUtils.h" +#include "Framework/CallbackService.h" +#include "Framework/ControlService.h" +#include "Framework/Task.h" +#include "Framework/CompletionPolicy.h" +#include "Framework/CompletionPolicyHelpers.h" +#include "Framework/DispatchPolicy.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/Variant.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/NameConf.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "Framework/CallbacksPolicy.h" + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// customize the completion policy +void customize(std::vector& policies) +{ + using o2::framework::CompletionPolicy; + using o2::framework::CompletionPolicyHelpers; + policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName("clusters-hmpid-write", CompletionPolicy::CompletionOp::Consume)); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + std::string keyvaluehelp("Semicolon separated key=value strings ..."); + workflowOptions.push_back(o2::framework::ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {keyvaluehelp}}); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +#include "Framework/runDataProcessing.h" +#include "HMPIDWorkflow/DigitsToClustersSpec.h" + +using namespace o2; +using namespace o2::framework; + +WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) +{ + WorkflowSpec specs; + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + DataProcessorSpec consumer = o2::hmpid::getDigitsToClustersSpec(); + specs.push_back(consumer); + return specs; +} diff --git a/Detectors/HMPID/workflow/src/digits-to-raw-stream-workflow.cxx b/Detectors/HMPID/workflow/src/digits-to-raw-stream-workflow.cxx index 1b619ab523e61..a200766162454 100644 --- a/Detectors/HMPID/workflow/src/digits-to-raw-stream-workflow.cxx +++ b/Detectors/HMPID/workflow/src/digits-to-raw-stream-workflow.cxx @@ -43,7 +43,7 @@ void customize(std::vector& workflowOptions) } #include "Framework/runDataProcessing.h" -#include "../include/HMPIDWorkflow/WriteRawFileSpec.h" +#include "HMPIDWorkflow/WriteRawFileSpec.h" using namespace o2; using namespace o2::framework; diff --git a/Detectors/HMPID/workflow/src/digits-to-root-workflow.cxx b/Detectors/HMPID/workflow/src/digits-to-root-workflow.cxx new file mode 100644 index 0000000000000..ce2987e67efb5 --- /dev/null +++ b/Detectors/HMPID/workflow/src/digits-to-root-workflow.cxx @@ -0,0 +1,71 @@ +// 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. + +/// \file digits-to-root-workflow.cxx +/// \author Antonio Franco - INFN Bari +/// \version 1.0 +/// \date 22 nov 2021 +/// + +#include "Framework/WorkflowSpec.h" +#include "Framework/DataSpecUtils.h" +#include "Framework/CallbackService.h" +#include "Framework/ControlService.h" +#include "Framework/Task.h" +#include "Framework/CompletionPolicy.h" +#include "Framework/CompletionPolicyHelpers.h" +#include "Framework/DispatchPolicy.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/Variant.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/NameConf.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "Framework/CallbacksPolicy.h" + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// customize the completion policy +void customize(std::vector& policies) +{ + using o2::framework::CompletionPolicy; + using o2::framework::CompletionPolicyHelpers; + policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName("digits-hmpid-root", CompletionPolicy::CompletionOp::Consume)); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + std::string keyvaluehelp("Semicolon separated key=value strings ..."); + workflowOptions.push_back(o2::framework::ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {keyvaluehelp}}); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); +} + +#include "Framework/runDataProcessing.h" +#include "HMPIDWorkflow/DigitsToRootSpec.h" + +using namespace o2; +using namespace o2::framework; + +WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) +{ + WorkflowSpec specs; + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + DataProcessorSpec consumer = o2::hmpid::getDigitsToRootSpec(); + specs.push_back(consumer); + + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + + return specs; +} diff --git a/Detectors/HMPID/workflow/src/raw-to-digits-workflow.cxx b/Detectors/HMPID/workflow/src/raw-to-digits-workflow.cxx index 429ac849cc63b..46e9c60302ce2 100644 --- a/Detectors/HMPID/workflow/src/raw-to-digits-workflow.cxx +++ b/Detectors/HMPID/workflow/src/raw-to-digits-workflow.cxx @@ -43,7 +43,7 @@ void customize(std::vector& workflowOptions) } #include "Framework/runDataProcessing.h" -#include "../include/HMPIDWorkflow/RawToDigitsSpec.h" +#include "HMPIDWorkflow/RawToDigitsSpec.h" using namespace o2; using namespace o2::framework;