diff --git a/core/data/tof/CbmTofHit.h b/core/data/tof/CbmTofHit.h
index bcb32c9b609aef17fbc47211b2335963fc84ca1a..d74dcc4277b89c0e0b43fa54fdc33db7c37d470f 100644
--- a/core/data/tof/CbmTofHit.h
+++ b/core/data/tof/CbmTofHit.h
@@ -82,6 +82,12 @@ public:
double GetCosPhi() const { return GetX() / GetRt(); }
double GetSinPhi() const { return GetY() / GetRt(); }
+ inline double Dist3D(CbmTofHit* pHit)
+ {
+ return sqrt((GetX() - pHit->GetX()) * (GetX() - pHit->GetX()) + (GetY() - pHit->GetY()) * (GetY() - pHit->GetY())
+ + (GetZ() - pHit->GetZ()) * (GetZ() - pHit->GetZ()));
+ }
+
/** Modifiers **/
void SetFlag(int32_t flag) { fFlag = flag; };
diff --git a/reco/detectors/tof/CbmTofCalibrator.cxx b/reco/detectors/tof/CbmTofCalibrator.cxx
index e38c4d3a677ebee0f626c39652bf0e1beeef1d65..6ffed6831f2bfb060209af7235c70ab775336705 100644
--- a/reco/detectors/tof/CbmTofCalibrator.cxx
+++ b/reco/detectors/tof/CbmTofCalibrator.cxx
@@ -16,8 +16,7 @@
#include "CbmTofAddress.h" // in cbmdata/tof
#include "CbmTofCell.h" // in tof/TofData
#include "CbmTofClusterizersDef.h"
-#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
-#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
+#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
#include "CbmTofDigiBdfPar.h" // in tof/TofParam
#include "CbmTofDigiPar.h" // in tof/TofParam
#include "CbmTofGeoHandler.h" // in tof/TofTools
@@ -30,6 +29,7 @@
#include <Logger.h>
// ROOT Classes and includes
+#include "TF1.h"
#include <TClonesArray.h>
#include <TDirectory.h>
#include <TFile.h>
@@ -37,10 +37,18 @@
#include <TGeoManager.h>
#include <TH1.h>
#include <TH2.h>
+#include <TH3.h>
#include <TMinuit.h>
#include <TProfile.h>
#include <TROOT.h>
+const int MaxShift = 10; // number of 50ps bins
+const double cLight = 29.9792; // in cm/ns
+const TString cBadChannelFile = "TofBadChannels.txt";
+const double dValidEdge = 1.0; // FIXME: constants in code
+
+std::vector<TH1*> hSvel;
+
CbmTofCalibrator::CbmTofCalibrator()
: fDigiMan(NULL)
, fTofClusterizer(NULL)
@@ -52,12 +60,21 @@ CbmTofCalibrator::CbmTofCalibrator()
, fhCalR0(NULL)
, fhCalDX0(NULL)
, fhCalDY0(NULL)
+ , fhCalCounterDt(NULL)
+ , fhCalCounterDy(NULL)
+ , fhCalChannelDt(NULL)
+ , fhCalChannelDy(NULL)
, fhCalTot()
, fhCalPosition()
, fhCalPos()
, fhCalTOff()
, fhCalTofOff()
+ , fhCalDelPos()
+ , fhCalDelTOff()
+ , fhCalCluTrms()
+ , fhCalCluSize()
, fhCalWalk()
+ , fhCalDtWalk()
, fhCorPos()
, fhCorTOff()
, fhCorTot()
@@ -128,8 +145,16 @@ InitStatus CbmTofCalibrator::Init()
return kFATAL;
}
+ fTofId = new CbmTofDetectorId_v21a();
+
fTrackletTools = new CbmTofTrackletTools(); // initialize tools
+ TString shcmd = Form("rm %s ", cBadChannelFile.Data());
+ gSystem->Exec(shcmd.Data());
+
+ shcmd = Form("touch %s ", cBadChannelFile.Data());
+ gSystem->Exec(shcmd.Data());
+
LOG(info) << "TofCalibrator initialized successfully";
return kSUCCESS;
}
@@ -164,12 +189,28 @@ Bool_t CbmTofCalibrator::CreateCalHist()
fhCalDX0 = new TH1D("hCalDX0", "Tracklet distance to nominal vertex; #DeltaX_0 [cm]", 100, -1.5, 1.5);
fhCalDY0 = new TH1D("hCalDY0", "Tracklet distance to nominal vertex; #DeltaY_0 [cm]", 100, -1.5, 1.5);
+ fhCalCounterDt = new TH1D("hCalCounterDt", "CalibCounterDt ; #Deltat [ns]", 100, -0.2, 0.2);
+ fhCalCounterDy = new TH1D("hCalCounterDy", "CalibCounterDy ; #Deltay [cm]", 100, -0.8, 0.8);
+ fhCalChannelDt = new TH1D("hCalChannelDt", "CalibChannelDt ; #Deltat [ns]", 100, -0.25, 0.25);
+ fhCalChannelDy = new TH1D("hCalChannelDy", "CalibChannelDy ; #Deltat [ns]", 100, -1.8, 1.8);
+
fhCalTot.resize(iNbDet);
fhCalPosition.resize(iNbDet);
fhCalPos.resize(iNbDet);
fhCalTOff.resize(iNbDet);
fhCalTofOff.resize(iNbDet);
+ fhCalDelPos.resize(iNbDet);
+ fhCalDelTOff.resize(iNbDet);
+ fhCalCluTrms.resize(iNbDet);
+ fhCalCluSize.resize(iNbDet);
fhCalWalk.resize(iNbDet);
+ fhCalDtWalk.resize(iNbDet);
+ fhCalWalkAv.resize(iNbDet);
+
+ fhCalXYTOff.resize(iNbDet);
+ fhCalXYTot.resize(iNbDet);
+ fhCalTotYWalk.resize(iNbDet);
+ fhCalTotYTOff.resize(iNbDet);
fDetIdIndexMap.clear();
@@ -186,6 +227,7 @@ Bool_t CbmTofCalibrator::CreateCalHist()
LOG(warning) << "No DigiPar for Det " << Form("0x%08x", iUCellId);
continue;
}
+ Double_t YMAX = TMath::Max(2., fChannelInfo->GetSizey()) * 0.75;
Double_t TotMax = 20.;
fhCalTot[iDetIndx] = new TH2F(
@@ -193,13 +235,18 @@ Bool_t CbmTofCalibrator::CreateCalHist()
Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; Strip []; Tot [a.u.]", iRpcId, iSmId, iSmType),
2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, 0., TotMax);
- Double_t YDMAX = TMath::Max(2., fChannelInfo->GetSizey()) * 1.;
+ fhCalXYTot[iDetIndx] =
+ new TH3F(Form("calXY_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId),
+ Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; Strip []; y [cm]; Tot [a.u.]", iRpcId, iSmId, iSmType),
+ 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 20, -YMAX,
+ YMAX, 100, 0., TotMax);
+
fhCalPosition[iDetIndx] =
new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_Position", iSmType, iSmId, iRpcId),
Form("Clu position of Rpc #%03d in Sm %03d of type %d; Strip []; ypos [cm]", iRpcId, iSmId, iSmType),
- fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YDMAX, YDMAX);
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YMAX, YMAX);
- YDMAX = 5;
+ Double_t YDMAX = 5;
fhCalPos[iDetIndx] =
new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_Pos", iSmType, iSmId, iRpcId),
Form("Clu pos deviation of Rpc #%03d in Sm %03d of type %d; "
@@ -209,31 +256,88 @@ Bool_t CbmTofCalibrator::CreateCalHist()
Double_t TSumMax = 2.;
//if(iSmType == 5) TSumMax *= 2.; // enlarge for diamond / beamcounter
- fhCalTOff[iDetIndx] = new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
- Form("Clu T0 deviation of Rpc #%03d in Sm %03d of type %d; "
- "Strip []; TOff [ns]",
- iRpcId, iSmId, iSmType),
- fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId),
- 99, -TSumMax, TSumMax);
- TSumMax = 50.;
+ fhCalTOff[iDetIndx] = new TH2F(
+ Form("cal_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
+ Form("Clu T0 deviation of Rpc #%03d in Sm %03d of type %d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 200, -TSumMax, TSumMax);
+
+ fhCalXYTOff[iDetIndx] = new TH3F(
+ Form("calXY_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
+ Form("XY T0 deviation of Rpc #%03d in Sm %03d of type %d; Strip []; y[cm]; TOff [ns]", iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YMAX, YMAX, 100,
+ -TSumMax * 0.5, TSumMax * 0.5);
+
+ TSumMax = 20.;
fhCalTofOff[iDetIndx] = new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_TofOff", iSmType, iSmId, iRpcId),
Form("Clu T0 deviation of Rpc #%03d in Sm %03d of type %d; "
"Strip []; TofOff [ns]",
iRpcId, iSmId, iSmType),
fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0,
- fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -TSumMax, TSumMax);
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 199, -TSumMax, TSumMax);
+
+ fhCalDelPos[iDetIndx] =
+ new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_DelPos", iSmType, iSmId, iRpcId),
+ Form("Clu position difference of Rpc #%03d in Sm %03d of type "
+ "%d; Strip []; #Deltaypos(clu) [cm]",
+ iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -10., 10.);
+
+ fhCalDelTOff[iDetIndx] =
+ new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_DelTOff", iSmType, iSmId, iRpcId),
+ Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
+ "[]; #DeltaTOff(clu) [ns]",
+ iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -0.6, 0.6);
+
+ fhCalCluTrms[iDetIndx] =
+ new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Trms", iSmType, iSmId, iRpcId),
+ Form("Clu Time RMS of Rpc #%03d in Sm %03d of type %d; Strip []; Trms [ns]", iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, 0., 0.5);
+
+ fhCalCluSize[iDetIndx] =
+ new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Size", iSmType, iSmId, iRpcId),
+ Form("Clu size of Rpc #%03d in Sm %03d of type %d; Strip []; size [strips]", iRpcId, iSmId, iSmType),
+ fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 16, 0.5, 16.5);
TSumMax = 0.8;
+ fhCalWalkAv[iDetIndx] =
+ new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_WalkAv", iSmType, iSmId, iRpcId),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_WalkAv; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId),
+ nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
+
fhCalWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
+ fhCalDtWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
+ fhCalTotYWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
+ fhCalTotYTOff[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
+
for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpcId); iCh++) {
fhCalWalk[iDetIndx][iCh].resize(2);
+ fhCalDtWalk[iDetIndx][iCh].resize(2);
+ fhCalTotYWalk[iDetIndx][iCh].resize(2);
+ fhCalTotYTOff[iDetIndx][iCh].resize(2);
for (Int_t iSide = 0; iSide < 2; iSide++) {
fhCalWalk[iDetIndx][iCh][iSide] =
new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
- Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot "
- "[a.u.]; #DeltaT [ns]",
- iSmType, iSmId, iRpcId, iCh, iSide),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId,
+ iRpcId, iCh, iSide),
+ nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
+ fhCalDtWalk[iDetIndx][iCh][iSide] =
+ new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_DtWalk", iSmType, iSmId, iRpcId, iCh, iSide),
+ Form("SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_DtWalk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId,
+ iCh, iSide),
nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
+
+ fhCalTotYWalk[iDetIndx][iCh][iSide] =
+ new TH3D(Form("calTotY_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
+ Form("YWalk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; y [cm]; #DeltaT [ns]", iSmType,
+ iSmId, iRpcId, iCh, iSide),
+ nbClWalkBinX, 0., TotMax, 20, -YMAX, YMAX, nbClWalkBinY, -0.4, 0.4);
+
+ fhCalTotYTOff[iDetIndx][iCh][iSide] =
+ new TH3D(Form("calTotY_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_TOff", iSmType, iSmId, iRpcId, iCh, iSide),
+ Form("YTot in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_TOff; Tot [a.u.]; y [cm]; #DeltaT [ns]", iSmType,
+ iSmId, iRpcId, iCh, iSide),
+ nbClWalkBinX, 0., TotMax, 20, -YMAX, YMAX, nbClWalkBinY, -TSumMax * 0.5, TSumMax * 0.5);
}
}
}
@@ -242,11 +346,13 @@ Bool_t CbmTofCalibrator::CreateCalHist()
}
static Int_t NevtH = 0;
-void CbmTofCalibrator::FillCalHist(CbmTofHit* pRef, Int_t /*iOpt*/, CbmEvent* /*tEvent*/)
+void CbmTofCalibrator::FillCalHist(CbmTofHit* pRef, Int_t iOpt, CbmEvent* tEvent)
{
- if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis!";
+ if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! for Event " << tEvent;
NevtH++;
+ if (NevtH == 1) LOG(info) << "FillCalHist: look for beam counter at " << pRef;
+
double dTAv = 0.;
if (pRef != NULL) dTAv = pRef->GetTime();
else {
@@ -347,18 +453,20 @@ void CbmTofCalibrator::FillCalHist(CbmTofHit* pRef, Int_t /*iOpt*/, CbmEvent* /*
LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
continue;
}
- gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
Double_t hitpos[3];
hitpos[0] = pHit->GetX();
hitpos[1] = pHit->GetY();
hitpos[2] = pHit->GetZ();
Double_t hlocal_p[3];
//TGeoNode* cNode=
- gGeoManager->GetCurrentNode();
- gGeoManager->MasterToLocal(hitpos, hlocal_p);
+ //gGeoManager->GetCurrentNode();
+ //gGeoManager->MasterToLocal(hitpos, hlocal_p);
+ fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
+
/*
- LOG(info)<<Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f ",
- NevtH,iSmType,iSm,iRpc,iCh,hlocal_p[1],pHit->GetTime() - pRef->GetTime());
+ LOG(info)<<Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f ",
+ NevtH,iSmType,iSm,iRpc,iCh,hlocal_p[1],pHit->GetTime() - dTAv, pRef->GetTime());
*/
fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
/*
@@ -366,19 +474,293 @@ void CbmTofCalibrator::FillCalHist(CbmTofHit* pRef, Int_t /*iOpt*/, CbmEvent* /*
fhCalTofOff[iDetIndx]->Fill((Double_t) iCh,
pHit->GetTime() - pRef->GetTime()); // ignore any dependence
*/
- fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ if (iOpt < 100) {
+ fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
+ fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ }
+ else {
+ double dTexp = pHit->GetR() * fTofFindTracks->GetTtTarg();
+ fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
+ fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
+ }
}
return;
}
+void CbmTofCalibrator::FillHitCalHist(CbmTofHit* pRef, const Int_t iOpt, CbmEvent* tEvent, TClonesArray* tTofHitsColl)
+{
+ if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! for Event " << tEvent;
+ if (NULL == tTofHitsColl) LOG(fatal) << "No access to TofHits! for Event " << tEvent;
+
+ NevtH++;
+ if (NevtH == 1)
+ LOG(warn) << "HitCalHist at " << NevtH << ", Opt = " << iOpt << ", pRef " << pRef << ", "
+ << Form("0x%08x", fTofClusterizer->GetBeamAddr());
+
+ double dTAv = 0.;
+ if (pRef != NULL) {
+ dTAv = pRef->GetTime();
+ //LOG(debug) << "HitCalHist got Tbeam " << dTAv << ", Opt = " << iOpt;
+ }
+ else {
+ if (-1 < fTofClusterizer->GetBeamAddr()) { // look for beam counter
+ if (NevtH == 1) LOG(warn) << Form("FillHitCalHist: look for beam counter 0x%08x", fTofClusterizer->GetBeamAddr());
+ Int_t iHit = 0;
+ //for (; iHit < (int) fTofFindTracks->GetNbHits(); iHit++) {
+ for (; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
+ //CbmTofHit* pHit = fTofFindTracks->GetHitPointer(iHit);
+ Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
+ CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
+ if (NULL == pHit) continue;
+ if ((pHit->GetAddress() & DetMask) == (fTofClusterizer->GetBeamAddr() & DetMask)) {
+ if (pRef != NULL) {
+ if (pHit->GetTime() < dTAv) pRef = pHit;
+ }
+ else {
+ pRef = pHit;
+ }
+ dTAv = pRef->GetTime();
+ if (NevtH == 1) LOG(warn) << "FillHitCalHist: found beam counter";
+ break;
+ }
+ }
+ if (iHit == tEvent->GetNofData(ECbmDataType::kTofHit)) return; // beam counter not present
+ }
+ else {
+ double dNAv = 0.;
+ for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
+ Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
+ CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
+ if (NULL == pHit) continue;
+ dTAv += pHit->GetTime();
+ dNAv++;
+ }
+ if (dNAv == 0) return;
+ dTAv /= dNAv;
+ }
+ }
+
+ for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
+ Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
+ CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
+ if (NULL == pHit) continue;
+
+ Int_t iDetId = (pHit->GetAddress() & DetMask);
+ Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
+ Int_t iSm = CbmTofAddress::GetSmId(iDetId);
+ Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
+
+ std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
+ if (it == fDetIdIndexMap.end()) {
+ LOG(info) << "detector TSR " << iSmType << iSm << iRpc << " not found in detector map";
+ continue; // continue for invalid detector index
+ }
+ Int_t iDetIndx = it->second; //fDetIdIndexMap[iDetId];
+
+ Int_t iChId = pHit->GetAddress();
+ Int_t iCh = CbmTofAddress::GetChannelId(iChId);
+ double dTexp = 0.;
+
+ /*
+ CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
+ if (NULL == fChannelInfo) {
+ LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
+ continue;
+ }
+ */
+
+ Double_t hitpos[3];
+ hitpos[0] = pHit->GetX();
+ hitpos[1] = pHit->GetY();
+ hitpos[2] = pHit->GetZ();
+ Double_t hlocal_p[3];
+
+ fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
+
+ LOG(debug) << Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f ", NevtH, iSmType, iSm, iRpc, iCh,
+ hlocal_p[1], pHit->GetTime() - dTAv, pRef->GetTime());
+
+ fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
+
+ if (iOpt < 100) {
+ LOG(debug) << "Fill " << fhCalTofOff[iDetIndx]->GetName();
+ if (pHit != pRef) {
+ fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
+ fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
+ }
+ }
+ else {
+ if (iSmType != 5) {
+ if (NULL != pRef) {
+ if ((pHit->GetAddress() & DetMask) != (pRef->GetAddress() & DetMask)) {
+ dTexp = pHit->Dist3D(pRef) / cLight; // Edge expected for speed of light,
+ if (pHit->GetZ() < pRef->GetZ()) dTexp *= -1;
+
+ LOG(debug) << Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f, Texp %12.1f ", NevtH, iSmType,
+ iSm, iRpc, iCh, hlocal_p[1], pHit->GetTime() - dTAv, pRef->GetTime(), dTexp);
+ }
+ }
+ else {
+ dTexp = pHit->GetR() / cLight; // Edge expected for speed of light, fTofFindTracks->GetTtTarg();
+ }
+ /*
+ LOG(debug)<<"Fill " << fhCalTofOff[iDetIndx]->GetName()<<" " << iHit <<": " << iCh << ", "
+ << pHit->GetTime() - dTAv <<" - " << dTexp;
+ */
+ if (dTexp != 0.) {
+ fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
+ fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv - dTexp);
+ fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
+ }
+ }
+ }
+
+ // fill CalDigi Tots and internal walks
+ if (dTexp != 0.) {
+ CbmMatch* digiMatch = fTofClusterizer->GetMatchIndexPointer(iHitIndex);
+ if (NULL == digiMatch) LOG(fatal) << "No digiMatch for Hit " << iHit << ", TSR " << iSmType << iSm << iRpc;
+ if (digiMatch->GetNofLinks() < fTofClusterizer->GetCluSizeMin()) continue;
+ double dMeanTimeSquared = 0.;
+ double dTotSum = 0.;
+ double dYSum = 0.;
+ for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // first loop over digis
+ CbmLink L0 = digiMatch->GetLink(iLink);
+ UInt_t iDigInd0 = L0.GetIndex();
+ UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
+ //LOG(debug1)<<" " << iDigInd0<<", "<<iDigInd1;
+ if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
+ const CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
+ const CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
+ if ((Int_t) pDig0->GetType() != iSmType) {
+ LOG(error) << Form(" Wrong Digi SmType %d - %d for Tofhit %d in iDetIndx "
+ "%d, Ch %d with %d strips at Indx %d, %d",
+ iSmType, (Int_t) pDig0->GetType(), iHit, iDetIndx, iCh,
+ fDigiBdfPar->GetNbChan(iSmType, iRpc), iDigInd0, iDigInd1);
+
+ for (Int_t iL = 0; iL < digiMatch->GetNofLinks(); iL++) { // loop over digis
+ int idx = (digiMatch->GetLink(iL)).GetIndex();
+ const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
+ LOG(info) << Form(" Digi %d, idx %d, TSRCS %d%d%d%02d%d ", iL, idx, (int) pDx->GetType(),
+ (int) pDx->GetSm(), (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
+ }
+
+ for (Int_t idx = 0; idx < (int) fTofCalDigiVec->size(); idx++) { // loop over digis
+ const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
+ LOG(info) << Form(" AllDigi %d, TSRCS %d%d%d%02d%d ", idx, (int) pDx->GetType(), (int) pDx->GetSm(),
+ (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
+ }
+
+ LOG(fatal) << " fhRpcCluTot: Digi 0 " << iDigInd0 << ": Ch " << pDig0->GetChannel() << ", Side "
+ << pDig0->GetSide() << ", StripSide " << (Double_t) iCh * 2. + pDig0->GetSide() << " Digi 1 "
+ << iDigInd1 << ": Ch " << pDig1->GetChannel() << ", Side " << pDig1->GetSide() << ", StripSide "
+ << (Double_t) iCh * 2. + pDig1->GetSide() << ", Tot0 " << pDig0->GetTot() << ", Tot1 "
+ << pDig1->GetTot();
+ } // end of type mismatch condition
+ fhCalTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
+ fhCalTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
+ if (digiMatch->GetNofLinks() > 2) {
+ double dTotStrip = pDig0->GetTot() + pDig1->GetTot();
+ dTotSum += dTotStrip;
+ double dYdigi = 0.5 * (pDig0->GetTime() - pDig1->GetTime()) * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ if (pDig0->GetSide() == 1) {
+ dYSum += (pDig0->GetTime() - pDig1->GetTime()) * dTotStrip; // weighted mean
+ }
+ else {
+ dYSum += (pDig1->GetTime() - pDig0->GetTime()) * dTotStrip;
+ dYdigi *= -1;
+ }
+ fhCalXYTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), dYdigi, pDig0->GetTot());
+ fhCalXYTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), dYdigi, pDig1->GetTot());
+ }
+ }
+ } // first loop over digis
+ if (digiMatch->GetNofLinks() > 2) {
+ double dYAv = dYSum / dTotSum;
+ double dYLoc = 0.;
+ for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // second loop over digis
+ CbmLink L0 = digiMatch->GetLink(iLink);
+ UInt_t iDigInd0 = L0.GetIndex();
+ UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
+ if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
+ const CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
+ const CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
+ double dTotStrip = pDig0->GetTot() + pDig1->GetTot();
+ double dTotWeight = 1 - dTotStrip / dTotSum;
+ //dTotWeight *= 0.5;
+ dTotWeight = 1.;
+ double dDeltaT = 0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->GetTime();
+ //if(NULL != fTofFindTracks) dDeltaT += fTofFindTracks->GetTOff(iDetId); //obsolete ??
+ fhCalDelTOff[iDetIndx]->Fill(pDig0->GetChannel(), dDeltaT * dTotWeight);
+ // TBD: position deviation in counter reference frame
+ Double_t dDelPos = 0.5 * (pDig0->GetTime() - pDig1->GetTime()) * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ if (0 == pDig0->GetSide()) dDelPos *= -1.;
+ fhCalDelPos[iDetIndx]->Fill(pDig0->GetChannel(), (dDelPos - fTofClusterizer->GetLocalY(pHit)));
+
+ int iCh0 = pDig0->GetChannel();
+ int iSide0 = pDig0->GetSide();
+ int iCh1 = pDig1->GetChannel();
+ int iSide1 = pDig1->GetSide();
+ if (iCh0 != iCh1) LOG(fatal) << "Inconsistent TofHit - FixIt";
+ fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(), dDeltaT * dTotWeight);
+ fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(), dDeltaT * dTotWeight);
+ dMeanTimeSquared += dDeltaT * dDeltaT;
+ if (pDig0->GetSide() == 1) { dYLoc = pDig0->GetTime() - pDig1->GetTime(); }
+ else {
+ dYLoc = pDig1->GetTime() - pDig0->GetTime();
+ }
+ dYLoc *= 0.5;
+ fhCalTotYWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(),
+ dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), dDeltaT);
+ fhCalTotYWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(),
+ dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), dDeltaT);
+
+ fhCalTotYTOff[iDetIndx][iCh0][iSide0]->Fill(
+ pDig0->GetTot(), dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), pHit->GetTime() - dTAv - dTexp);
+ fhCalTotYTOff[iDetIndx][iCh1][iSide1]->Fill(
+ pDig1->GetTot(), dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), pHit->GetTime() - dTAv - dTexp);
+
+ Double_t dDeltaDT = dYLoc - dYAv;
+ fhCalPos[iDetIndx]->Fill(iCh0, dDeltaDT * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
+
+ fhCalDtWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(), (2. * pDig0->GetSide() - 1) * dDeltaDT);
+ fhCalDtWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(), (2. * pDig1->GetSide() - 1) * dDeltaDT);
+ }
+ } // second loop over digis
+ } // if (digiMatch->GetNofLinks()>2)
+
+ fhCalCluSize[iDetIndx]->Fill((Double_t) iCh, digiMatch->GetNofLinks() / 2);
+ if (digiMatch->GetNofLinks() > 2) {
+ double dVar = dMeanTimeSquared / (digiMatch->GetNofLinks() / 2 - 1);
+ double dTrms = TMath::Sqrt(dVar);
+ fhCalCluTrms[iDetIndx]->Fill((Double_t) iCh, dTrms);
+ }
+ } // iOpt<100
+ }
+ return;
+} // -- FillHitCalHist
+
+
static Int_t NevtT = 0;
void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* tEvent)
{
NevtT++;
- if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis!";
+ if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! Opt " << iOpt << ", Event" << tEvent;
// fill deviation histograms on walk level
+ /*
+ LOG(info)<<"Calibrator " << NevtT <<": Tt " << pTrk->GetTt()
+ <<", R0: " << pTrk->GetR0() <<", bB "
+ << pTrk->ContainsAddr(CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 5));
+ */
if (pTrk->GetTt() < 0) return; // take tracks with positive velocity only
+ if (pTrk->GetNofHits() < fTofFindTracks->GetNReqStations()) return;
+
+ fhCalDX0->Fill(pTrk->GetFitX(0.));
+ fhCalDY0->Fill(pTrk->GetFitY(0.));
+
if (fbBeam && !pTrk->ContainsAddr(CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 5)))
return; // request beam counter hit for calibration
@@ -387,9 +769,8 @@ void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* t
fhCalR0->Fill(pTrk->GetR0());
if (pTrk->GetR0() > fdR0Lim) return;
}
- fhCalDX0->Fill(pTrk->GetFitX(0.));
- fhCalDY0->Fill(pTrk->GetFitY(0.));
- if (iOpt > 70) HstDoublets(pTrk); // Fill Doublets histograms
+
+ //if (iOpt > 70) HstDoublets(pTrk); // Fill Doublets histograms
//double dEvtStart=fEvtHeader->GetEventTime();
for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
@@ -410,24 +791,37 @@ void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* t
LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
continue;
}
- gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
Double_t hitpos[3];
hitpos[0] = pHit->GetX();
hitpos[1] = pHit->GetY();
hitpos[2] = pHit->GetZ();
Double_t hlocal_p[3];
//TGeoNode* cNode=
- gGeoManager->GetCurrentNode();
- gGeoManager->MasterToLocal(hitpos, hlocal_p);
+ //gGeoManager->GetCurrentNode();
+ //gGeoManager->MasterToLocal(hitpos, hlocal_p);
+ fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
+
hitpos[0] = pTrk->GetFitX(pHit->GetZ());
hitpos[1] = pTrk->GetFitY(pHit->GetZ());
Double_t hlocal_f[3];
- gGeoManager->MasterToLocal(hitpos, hlocal_f);
- fhCalPosition[iDetIndx]->Fill((Double_t) iCh,
- hlocal_p[1]); // transformed into LRF
- fhCalPos[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1] - hlocal_f[1]); // transformed into LRF
- fhCalTOff[iDetIndx]->Fill((Double_t) iCh,
- pHit->GetTime() - pTrk->GetFitT(pHit->GetZ())); // residuals transformed into LRF
+ //gGeoManager->MasterToLocal(hitpos, hlocal_f);
+ fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_f);
+ if (true) { //0==fDigiBdfPar->GetChanOrient(iSmType, iRpc)) { // default
+ fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
+ fhCalPos[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1] - hlocal_f[1]); // transformed into LRF
+ }
+ else { // Counters rotated with respect to module axis, strips along x-axis
+ fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[0]); // transformed into LRF
+ fhCalPos[iDetIndx]->Fill((Double_t) iCh, hlocal_p[0] - hlocal_f[0]); // transformed into LRF
+ }
+ //if(fTofFindTracks->GetTtTarg()==0.033) pTrk->SetTt(fTofFindTracks->GetTtTarg()); //for calibration mode 1,2,31
+ //double dTexp = pTrk->GetFitT(pHit->GetZ()); // residuals transformed into LRF
+ //double dTexp = fTrackletTools->GetTexpected(pTrk, pHit->GetAddress()&DetMask, pHit, fTofFindTracks->GetTtLight());
+ double dTexp = fTrackletTools->GetTexpected(pTrk, -1, pHit, fTofFindTracks->GetTtLight());
+ fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTexp);
+ fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1],
+ pHit->GetTime() - dTexp); // residuals transformed into LRF
//fhCalTOff[iDetIndx]->Fill((Double_t)iCh,fTrackletTools->GetTdif(pTrk, iDetId, pHit));
Int_t iEA = pTrk->GetTofHitIndex(iHit);
@@ -486,22 +880,25 @@ void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* t
continue;
}
if (iCh0 > (Int_t) fhCalWalk[iDetIndx].size()) {
- LOG(error) << "Invalid Ch0 " << iCh0;
+ LOG(error) << "Invalid Ch0 " << iCh0 << " for detIndx" << iDetIndx << ", size " << fhCalWalk[iDetIndx].size();
continue;
}
if (iSide0 > (Int_t) fhCalWalk[iDetIndx][iCh0].size()) {
- LOG(error) << "Invalid Side0 " << iSide0;
+ LOG(error) << "Invalid Side0 " << iSide0 << " for "
+ << " for detIndx" << iDetIndx << ", size " << fhCalWalk[iDetIndx][iCh0].size();
continue;
}
Int_t iCh1 = tDigi1->GetChannel();
Int_t iSide1 = tDigi1->GetSide();
if (iCh1 > (Int_t) fhCalWalk[iDetIndx].size()) {
- LOG(error) << "Invalid Ch1 " << iCh1;
+ LOG(error) << "Invalid Ch1 " << iCh1 << " for "
+ << " for detIndx" << iDetIndx;
continue;
}
if (iSide1 > (Int_t) fhCalWalk[iDetIndx][iCh1].size()) {
- LOG(error) << "Invalid Side1 " << iSide1;
+ LOG(error) << "Invalid Side1 " << iSide1 << " for "
+ << " for detIndx" << iDetIndx;
continue;
}
@@ -519,7 +916,7 @@ void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* t
LOG(warn) << "CMPY for TSRC " << iSmType << iSm << iRpc << iCh0 << ": " << hlocal_f[1] << ", " << hlocal_p[1]
<< ", " << hlocal_d[1] << ", TOT: " << tDigi0->GetTot() << " " << tDigi1->GetTot();
}
- Int_t iWalkMode = (iOpt - iOpt % 100) / 100;
+ Int_t iWalkMode = 0; //(iOpt - iOpt % 100) / 100;
switch (iWalkMode) {
case 1:
fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(
@@ -549,14 +946,20 @@ void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* t
break;
case 0: {
- Double_t dDeltaT = 0.5 * (tDigi0->GetTime() + tDigi1->GetTime()) + fTofFindTracks->GetTOff(iDetId)
- - pTrk->GetFitT(pHit->GetZ());
+ Double_t dDeltaT =
+ 0.5 * (tDigi0->GetTime() + tDigi1->GetTime()) + fTofFindTracks->GetTOff(iDetId) - pHit->GetTime();
+
+ fhCalWalkAv[iDetIndx]->Fill(tDigi0->GetTot(), dDeltaT);
+ fhCalWalkAv[iDetIndx]->Fill(tDigi1->GetTot(), dDeltaT);
fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi0->GetTot(), dDeltaT);
fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi1->GetTot(), dDeltaT);
if (iSmType == 5 || iSmType == 8) { // symmetrize for Pad counters
fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi1->GetTot(), dDeltaT);
fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi0->GetTot(), dDeltaT);
}
+ Double_t dDeltaDT = (hlocal_d[1] - hlocal_p[1]) / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ fhCalDtWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi0->GetTot(), (2. * tDigi0->GetSide() - 1) * dDeltaDT);
+ fhCalDtWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi1->GetTot(), (2. * tDigi1->GetSide() - 1) * dDeltaDT);
} break;
}
}
@@ -630,126 +1033,779 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
switch (iOpt0) {
case 0: // none
- break;
- case 1: // update channel mean
{
- LOG(info) << "Update Offsets for TSR " << iSmType << iSm << iRpc;
+ double* dRes;
+ switch (iOpt1) {
+ case 10: {
+ LOG(info) << "Equalize velocities with const threshold counterwise for " << fhCorTOff[iDetIndx]->GetName()
+ << ", MeanTOff " << fhCorTOff[iDetIndx]->GetMean();
+ TString hname2 = Form("cal_SmT%d_sm%03d_rpc%03d_TOff", iSmType, iSm, iRpc);
+ //TH2* h2 = (TH2*) gROOT->FindObjectAny(hname2);
+ TH2* h2 = fhCalTOff[iDetIndx];
+ TH2* h2W = fhCalTofOff[iDetIndx];
+ if (NULL == h2) {
+ LOG(warn) << "Calibration data not available from " << hname2;
+ continue;
+ }
+
+ double dResIni[1] = {0.};
+ dRes = &dResIni[0];
+ Int_t NStrips = h2->GetNbinsX();
+ TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
+ TH1* h2y = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
+ double dEdge = 0.;
+ if (iSmType != 5) {
+ dRes = find_tofedge((const char*) (h2Wy->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ dEdge = dRes[0];
+ }
+ else {
+ dEdge = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
+ }
+
+ /*
+ LOG(info) << h2Wy->GetName() << ": Coarse Max at " << h2Wy->GetMaximumBin()
+ << ", " << dEdge << ", " << h2y->GetBinCenter(1)
+ << ", " << h2y->GetBinCenter(h2y->GetNbinsX());
+ */
+ if (dEdge < h2y->GetBinCenter(1) || dEdge > h2y->GetBinCenter(h2y->GetNbinsX())) {
+ dRes[0] = dEdge; // - fTofClusterizer->GetEdgeTbias();
+ LOG(info) << h2Wy->GetName() << ": Coarse TOff shift by " << dRes[0];
+ }
+ else {
+ if (iSmType != 5) {
+ dRes = find_tofedge((const char*) (h2y->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ if (TMath::Abs(dRes[0]) < 0.5) {
+ //LOG(info) << "Now fit edge above noise, dT " << dRes[0] <<", T "<< iSmType;
+ double dTfind = dRes[0];
+ double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
+ dRes = fit_tofedge((const char*) (h2y->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
+ if (isnan(dRes[0])) dRes[0] = dTfind;
+ }
+ }
+ }
+ //LOG(info) << h2y->GetName() << ": shift TOff by " << dRes[0];
+ fhCalChannelDt->Fill(dRes[0]);
+ for (Int_t i = 0; i < NStrips; i++) {
+ fhCorTOff[iDetIndx]->SetBinContent(i + 1, fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0]);
+ }
+ } break;
+
+ case 15:
+ case 12:
+ case 11: {
+ LOG(info) << "Equalize velocities with chi2 deviations for " << fhCorTOff[iDetIndx]->GetName();
+ TString hname2 = Form("cal_SmT%d_sm%03d_rpc%03d_TOff", iSmType, iSm, iRpc);
+
+ double dTRefSum = 0.;
+ double dRefNormSum = 0.;
+ double dNorm = 0.;
+
+ TH2* h2W = fhCalTofOff[iDetIndx];
+
+ if (h2W->GetEntries() < MINCTS) continue;
+
+ TH2* h2 = fhCalTOff[iDetIndx];
+ Int_t NStrips = h2->GetNbinsX();
+ TH1* h2yAv;
+ double dResIni[1] = {0.};
+ dRes = &dResIni[0];
+
+ TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
+ h2yAv = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
+ double dEdge = 0.;
+
+ double dMean = TruncatedMeanY(h2W); //h2Wy->GetMean();
+ if (iSmType == 5) {
+ dEdge = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
+ dNorm = h2Wy->GetBinContent(h2Wy->GetMaximumBin());
+ }
+ else {
+ dRes = find_tofedge((const char*) (h2Wy->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ dEdge = dRes[0];
+ if ((dMean - dEdge) > 10.) { //FIXME: constant in code
+ LOG(warn) << h2Wy->GetName() << ": wrong edge detected " << dEdge << ", Mean " << dMean;
+ dEdge = dMean;
+ }
+ }
+ /*
+ LOG(info) << h2Wy->GetName() << ": Coarse shift at " << h2Wy->GetMaximumBin()
+ << ", Edge " << dEdge
+ << ", " << h2yAv->GetBinCenter(1)
+ << ", " << h2yAv->GetBinCenter(h2yAv->GetNbinsX());
+ */
+ if (dEdge < h2yAv->GetBinCenter(1) || dEdge > h2yAv->GetBinCenter(h2yAv->GetNbinsX())) {
+ dRes[0] = dEdge; // - fTofClusterizer->GetEdgeTbias();
+ //LOG(info) << h2Wy->GetName() << ": Coarse TOff shift by " << dRes[0];
+ }
+ else {
+ if (iSmType != 5) {
+ dRes = find_tofedge((const char*) (h2yAv->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ if (TMath::Abs(dRes[0]) < dValidEdge) {
+ double dTfind = dRes[0];
+ double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
+ dRes = fit_tofedge((const char*) (h2yAv->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
+ if (isnan(dRes[0])) dRes[0] = dTfind;
+ LOG(info) << h2yAv->GetName() << ": fitted AvEdge up to " << dTmax << ", Res " << dRes[0]
+ << ", FindE " << dTfind;
+ }
+ }
+ }
+
+ if (iSmType != 5) {
+ fhCalCounterDt->Fill(dRes[0] - fTofClusterizer->GetEdgeTbias());
+ if (TMath::Abs(dRes[0] - fTofClusterizer->GetEdgeTbias()) > -fhCalCounterDt->GetBinCenter(1)) {
+ LOG(warn) << "Channel ooCR: " << fhCorTOff[iDetIndx]->GetName() << ", " << dRes[0] << ", "
+ << fTofClusterizer->GetEdgeTbias() << ", " << -fhCalCounterDt->GetBinCenter(1);
+ }
+ }
+ else {
+ fhCalCounterDt->Fill(dEdge);
+ }
+
+ const double dResAv = dRes[0];
+ const int nValues = MaxShift * 2 + 1;
+ const double dRangeFac = 1.;
+ double chi2Shift[NStrips][nValues];
+ double chi2Minimum[NStrips];
+ double chi2MinShift[NStrips];
+ double chi2LowMinNeighbor[NStrips];
+ double chi2LowMinNeighborShift[NStrips];
+ for (Int_t i = 0; i < NStrips; i++) {
+ double dTminShift = 0.;
+ double dEdgei = 0.;
+ double dMeani = 0.;
+ TH1* h2y = h2->ProjectionY(Form("%s_py%d", h2->GetName(), i), i + 1, i + 1);
+ if (h2y->GetEntries() < MINCTS) continue;
+ LOG(info) << h2y->GetName() << ": inspect with iOpt1 = " << iOpt1;
+ if (kTRUE) {
+ if (iOpt1 == 11 || iOpt1 == 15) {
+ //check peak position
+ TH1* h2Wyi = h2W->ProjectionY(Form("%s_py%d", h2W->GetName(), i), i + 1, i + 1);
+ if (iSmType == 5) {
+ dEdgei = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
+ dNorm = h2Wy->GetBinContent(h2Wy->GetMaximumBin());
+ }
+ else {
+ dRes = find_tofedge((const char*) (h2Wyi->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ dEdgei = dRes[0];
+ }
+
+ if (dEdgei < dRangeFac * h2y->GetBinCenter(1)
+ || dEdgei > dRangeFac * h2y->GetBinCenter(h2y->GetNbinsX())) {
+ dMeani = h2Wyi->GetMean();
+ LOG(info) << h2Wyi->GetName() << " ooR: " << dRes[0] << ", E " << dEdgei << ", Av " << dResAv
+ << ", Meani " << dMeani << ", MeanY " << dMean;
+ if ((dMeani - dEdgei) > 10.) { //FIXME: constant in code
+ LOG(warn) << h2Wyi->GetName() << ": invalid edge detected " << dEdgei << ", Mean " << dMeani;
+ dRes[0] = dMeani - dMean;
+ }
+ //dRes[0]=dEdgei
+ }
+ else {
+ if (iSmType != 5) {
+ dRes = find_tofedge((const char*) (h2y->GetName()), fTofClusterizer->GetEdgeThr(),
+ fTofClusterizer->GetEdgeLen());
+ LOG(info) << h2y->GetName() << " CheckValidEdge: " << dRes[0] << ", " << dValidEdge;
+ if (TMath::Abs(dRes[0]) < dValidEdge) {
+ switch (iOpt1) {
+ case 11: {
+ double dTfind = dRes[0];
+ double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
+ dRes = fit_tofedge((const char*) (h2y->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
+ if (isnan(dRes[0])) dRes[0] = dTfind;
+ LOG(info) << h2y->GetName() << ": fit edge up to " << dTmax << ", Res " << dRes[0]
+ << ", FindE " << dTfind;
+ } break;
+ case 15:
+ chi2Minimum[i] = 1.E6;
+ chi2LowMinNeighbor[i] = 1.E6;
+ int idx = 0;
+ for (int j = -MaxShift; j < MaxShift + 1; j++) {
+ chi2Shift[i][idx] = CalcChi2(h2yAv, h2y, j);
+ if (chi2Shift[i][idx] < chi2Minimum[i]) {
+ chi2LowMinNeighbor[i] = chi2Minimum[i];
+ chi2LowMinNeighborShift[i] = chi2MinShift[i];
+ chi2Minimum[i] = chi2Shift[i][idx];
+ chi2MinShift[i] = j * h2y->GetBinWidth(1);
+ }
+ else {
+ if (chi2Shift[i][idx] < chi2LowMinNeighbor[i]) {
+ chi2LowMinNeighbor[i] = chi2Shift[i][idx];
+ chi2LowMinNeighborShift[i] = j * h2y->GetBinWidth(1);
+ }
+ }
+ idx++;
+ }
+ double chi2sum = chi2Minimum[i] + chi2LowMinNeighbor[i];
+ dTminShift = (chi2MinShift[i] * (chi2sum - chi2Minimum[i])
+ + chi2LowMinNeighborShift[i] * (chi2sum - chi2LowMinNeighbor[i]))
+ / chi2sum;
+ //dTminShift=chi2MinShift[i]; // take Min only, for debugging
+ LOG(info) << h2y->GetName() << ": Chi2 shift " << dTminShift << ", Res " << dRes[0]
+ << ", ResAv " << dResAv;
+ dRes[0] = dResAv;
+ break;
+ }
+ }
+ else {
+ LOG(info) << h2y->GetName() << ": stick to coarse offset " << dEdgei << ", Res " << dRes[0];
+ dRes[0] = dEdgei; // use offset from wide histo
+ }
+ }
+ else { // beam counter
+ dEdgei = h2y->GetBinCenter(h2y->GetMaximumBin());
+ dNorm = h2y->GetBinContent(h2y->GetMaximumBin());
+ }
+
+ //LOG(info) << h2Wyi->GetName() << " TminShift: " << dTminShift;
+ }
+ }
+ }
+ if (iSmType == 5) {
+ LOG(info) << "Update StartCounterCalib " << i << ": " << dRes[0] << ", " << dBeamTOff << ", "
+ << dTminShift << ", " << dEdgei << ", N " << dNorm;
+ dRes[0] = -dEdgei;
+ }
+ else {
+ dTminShift -= fTofClusterizer->GetEdgeTbias(); // shift to physical scale
+ }
+ fhCalChannelDt->Fill(dRes[0] + dTminShift);
+ if (TMath::Abs(dRes[0] + dTminShift) > -fhCalChannelDt->GetBinCenter(1)) {
+ LOG(warn) << "Channel ooHR: " << h2Wy->GetName() << i << ", " << dRes[0] << ", " << dTminShift;
+ }
+ double dVal = fhCorTOff[iDetIndx]->GetBinContent(i + 1);
+ if (isnan(dVal)) fhCorTOff[iDetIndx]->SetBinContent(i + 1, 0.); // fix NANs
+
+ LOG(info) << h2Wy->GetName() << i << " CorTOff: " << dRes[0] << ", " << dResAv << ", " << dTminShift
+ << ", " << fhCorTOff[iDetIndx]->GetBinContent(i + 1) << " -> "
+ << fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0] + dTminShift;
+
+ fhCorTOff[iDetIndx]->SetBinContent(i + 1,
+ (fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0] + dTminShift));
+ if (iSmType == 5) {
+ dTRefSum += fhCorTOff[iDetIndx]->GetBinContent(i + 1) * dNorm;
+ dRefNormSum += dNorm;
+ }
+ //LOG(info)<<fhCorTOff[iDetIndx]->GetName() << ": new TOff " << fhCorTOff[iDetIndx]->GetBinContent(i + 1);
+ } //for (Int_t i = 0; i < NStrips; i++)
+ if (iSmType == 5) { // avoid runaway
+ double dTRefAv = dTRefSum / dRefNormSum;
+ LOG(info) << "Shift average ref time to 0: " << dTRefAv;
+ for (Int_t i = 0; i < NStrips; i++) {
+ fhCorTOff[iDetIndx]->SetBinContent(i + 1, (fhCorTOff[iDetIndx]->GetBinContent(i + 1) - dTRefAv));
+ }
+ }
+ } break;
+
+ case 14:
+ case 13: {
+ LOG(info) << "TofCalibrator: skip TOff, calibrate hit positions, check YBox, Opt = " << iOpt;
+ } break;
+ case 16: {
+ LOG(info) << "TofCalibrator: determine counter time offsets with cosmics, Opt = " << iOpt;
+ TH2* h2 = fhCalTOff[iDetIndx];
+ TH2* h2W = fhCalTofOff[iDetIndx];
+ if (NULL == h2) {
+ LOG(warn) << "Calibration data not available for detx " << iDetIndx;
+ continue;
+ }
+ Int_t NStrips = h2->GetNbinsX();
+ TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
+ TH1* h2y = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
+
+ Double_t dCtsWy = h2Wy->GetEntries();
+ Double_t dDt = 0;
+ if (dCtsWy > MINCTS) {
+ double dFMean = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
+ double dFLim = 0.5; // CAUTION, fixed numeric value
+ double dBinSize = h2Wy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = h2Wy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1);
+ if (TMath::Abs(dDt) < 0.9 * h2y->GetXaxis()->GetXmax()) {
+ dFMean = h2y->GetBinCenter(h2y->GetMaximumBin());
+ dFLim = 0.5; // CAUTION, fixed numeric value
+ dBinSize = h2y->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ fRes = h2y->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1);
+ }
+ }
+ fhCalCounterDt->Fill(dDt);
+ }
+ LOG(info) << "Update cor hist " << fhCorTOff[iDetIndx]->GetName() << " by " << dDt << ", " << dBeamTOff;
+ for (int iCh = 0; iCh < fhCorTOff[iDetIndx]->GetNbinsX(); iCh++) {
+ fhCalChannelDt->Fill(dDt);
+ fhCorTOff[iDetIndx]->SetBinContent(iCh + 1, fhCorTOff[iDetIndx]->GetBinContent(iCh + 1) + dDt);
+ }
+ }
+ } break;
+
+ case 17: {
+ LOG(info) << "TofCalibrator: determine channel time offsets with cosmics, Opt = " << iOpt;
+ TH2* h2 = fhCalTOff[iDetIndx];
+ TH2* h2W = fhCalTofOff[iDetIndx];
+ if (NULL == h2) {
+ LOG(warn) << "Calibration data not available for detx " << iDetIndx;
+ continue;
+ }
+ Int_t NStrips = h2->GetNbinsX();
+ for (int i = 0; i < NStrips; i++) {
+ TH1* h2Wy = h2W->ProjectionY(Form("%s_py_%02d", h2W->GetName(), i), i + 1, i + 1);
+ TH1* h2y = h2->ProjectionY(Form("%s_py_%02d", h2->GetName(), i), i + 1, i + 1);
+
+ Double_t dCtsWy = h2Wy->GetEntries();
+ Double_t dDt = 0;
+ if (dCtsWy > MINCTS) {
+ double dFMean = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
+ double dFLim = 0.5; // CAUTION, fixed numeric value
+ double dBinSize = h2Wy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = h2Wy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ Int_t iFitStatus = fRes;
+ //LOG(warn)<<h2Wy->GetName() << ": "<< iFitStatus << " " <<fRes;
+ if (iFitStatus != -1)
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1);
+ //LOG(warn)<<h2Wy->GetName() << " Dt "<< dDt;
+ if (h2y->GetEntries() > MINCTS)
+ if (TMath::Abs(dDt) < 0.9 * h2y->GetXaxis()->GetXmax()) {
+ dFMean = h2y->GetBinCenter(h2y->GetMaximumBin());
+ dFLim = 0.5; // CAUTION, fixed numeric value
+ dBinSize = h2y->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ fRes = h2y->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1);
+ }
+ }
+ fhCalChannelDt->Fill(dDt);
+ fhCorTOff[iDetIndx]->SetBinContent(i + 1, fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dDt);
+ }
+ }
+ else {
+ LOG(warn) << "Too few counts in " << h2Wy->GetName() << ": " << dCtsWy;
+ }
+ }
+ } break;
+
+ case 18: {
+ LOG(info) << "TofCalibrator: determine channel walks with cosmics, Opt = " << iOpt << " for TSR "
+ << iSmType << iSm << iRpc;
+ //if (iSmType == 5) continue; // no walk correction for beam counter up to now
+ const Double_t MinCounts = 10.;
+ Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
+ TH1* hCY = fhCalPos[iDetIndx]->ProjectionY(Form("%s_py_%d", fhCalPos[iDetIndx]->GetName(), iCh),
+ iCh + 1, iCh + 1);
+ if (hCY->GetEntries() > 1)
+ fhCalChannelDy->Fill(hCY->GetRMS() / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
+ //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
+ TProfile* hpW0 = fhCalWalk[iDetIndx][iCh][0]->ProfileX(); // mean deviation
+ TH1* hCW0 = fhCalWalk[iDetIndx][iCh][0]->ProjectionX(); // contributing counts
+ TProfile* hpW1 = fhCalWalk[iDetIndx][iCh][1]->ProfileX(); // mean deviation
+ TH1* hCW1 = fhCalWalk[iDetIndx][iCh][1]->ProjectionX(); // contributing counts
+ if (hCW0->GetEntries() == 0 || hCW1->GetEntries() == 0) {
+ //LOG(info) << "No entries in " << hCW0->GetName();
+ continue;
+ }
+ Double_t dCorT = 0;
+ for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][0]->GetNbinsX(); iBin++) {
+ Double_t dCts0 = hCW0->GetBinContent(iBin + 1);
+ Double_t dCts1 = hCW1->GetBinContent(iBin + 1);
+ Double_t dWOff0 = fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin + 1); // current value
+ Double_t dWOff1 = fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin + 1); // current value
+ if (iBin > 0 && dCts0 == 0 && dCts1 == 0) {
+ fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1,
+ fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin));
+ fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1,
+ fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin));
+ }
+ dCorT = 0.;
+ if (dCts0 > MinCounts && dCts1 > MinCounts) {
+ dCorT = 0.5 * (hpW0->GetBinContent(iBin + 1) + hpW1->GetBinContent(iBin + 1));
+ fhCalChannelDt->Fill(dCorT);
+ }
+ fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1, dWOff0 + dCorT); //set new value
+ fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1, dWOff1 + dCorT); //set new value
+ }
+ }
+ } break;
+
+ default:
+ LOG(info) << "TofCalibrator: unknown option " << iOpt;
+ return kTRUE;
+ break;
+ } // switch iOpt1 end
+
+ // re-adjust positions
+ if (fhCalPosition[iDetIndx]->GetEntries() == 0) continue;
+ double dYShift = 0;
+ TH1* hCalP = fhCalPosition[iDetIndx]->ProjectionX();
+ TH1* hCalPos_py = fhCalPosition[iDetIndx]->ProjectionY(); // full counter
+ dYShift = hCalPos_py->GetMean();
+ int iChId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
+ CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
+ if (NULL == fChannelInfo) LOG(fatal) << Form("invalid ChannelInfo for 0x%08x", iChId);
+ double dLen = fChannelInfo->GetSizey() * fTofClusterizer->GetModifySigvel();
+ fTofClusterizer->fit_ybox(hCalPos_py, dLen);
+ TF1* ff = hCalPos_py->GetFunction("YBox"); // box fit for monitoring signal velocity
+ if (NULL != ff) {
+ LOG(info) << "FRes YBox " << hCalPos_py->GetEntries() << " entries in TSR " << iSmType << iSm << iRpc
+ << ", chi2 " << ff->GetChisquare() / ff->GetNDF()
+ << Form(", striplen (%5.2f): %7.2f +/- %5.2f, pos "
+ "res %5.2f +/- %5.2f at y_cen = %5.2f +/- %5.2f",
+ dLen * fTofClusterizer->GetModifySigvel(), 2. * ff->GetParameter(1),
+ 2. * ff->GetParError(1), ff->GetParameter(2), ff->GetParError(2), ff->GetParameter(3),
+ ff->GetParError(3));
+ //dYShift = ff->GetParameter(3); // update common shift by box fit
+ if (iOpt1 == 14) { // Update signal velocity
+ if (hSvel[iSmType] != NULL) {
+ int iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ double dSscal = hSvel[iSmType]->GetBinContent(iSm * iNbRpc + iRpc + 1);
+ double dRatio = dLen * fTofClusterizer->GetModifySigvel() / 2. / ff->GetParameter(1);
+ LOG(info) << "Modify SigVel modifier for TSR " << iSmType << iSm << iRpc << ": " << dSscal << ", "
+ << dRatio << " -> " << dSscal * dRatio;
+ dSscal = dSscal * dRatio;
+ if (dSscal < 0.8) dSscal = 0.8;
+ if (dSscal > 1.2) dSscal = 1.2;
+ hSvel[iSmType]->SetBinContent(iSm * iNbRpc + iRpc + 1,
+ dSscal); // does not work in input file directory
+ }
+ }
+ }
+
+ double dThr = hCalPos_py->GetBinContent(hCalPos_py->GetMaximumBin()) / 10;
+ if (dThr < 3.) dThr = 3.;
+ dRes = fTofClusterizer->find_yedges((const char*) (hCalPos_py->GetName()), dThr, fChannelInfo->GetSizey());
+ if (TMath::Abs(dRes[0] - fChannelInfo->GetSizey()) / fChannelInfo->GetSizey() < 0.1) { dYShift = dRes[1]; }
+
+ if (iOpt1 != 16 && iOpt1 != 12) {
+ for (Int_t iBin = 0; iBin < fhCorPos[iDetIndx]->GetNbinsX(); iBin++) {
+ Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
+ Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
+ if (dCtsP > MINCTS) {
+ TH1* hPos_py = fhCalPosition[iDetIndx]->ProjectionY(
+ Form("%s_py_%d", fhCalPosition[iDetIndx]->GetName(), iBin), iBin + 1, iBin + 1);
+ dYShift = hPos_py->GetMean();
+ dThr = hPos_py->GetBinContent(hPos_py->GetMaximumBin()) / 10;
+ if (dThr < 3.) dThr = 3.;
+ dRes = fTofClusterizer->find_yedges((const char*) (hPos_py->GetName()), dThr, fChannelInfo->GetSizey());
+
+ LOG(info) << Form("EdgeY for %s, TSR %d%d%d: DY %5.2f, Len %5.2f, Size %5.2f ", hPos_py->GetName(),
+ iSmType, iSm, iRpc, dRes[1], dRes[0],
+ fChannelInfo->GetSizey() * fTofClusterizer->GetModifySigvel());
+
+ if (TMath::Abs(dRes[0] - fChannelInfo->GetSizey()) / fChannelInfo->GetSizey() < 0.1) {
+ dYShift = dRes[1];
+ }
+ else {
+ dYShift = hPos_py->GetMean();
+ }
+ // apply correction
+ //LOG(info)<<Form("UpdateCalPos TSRC %d%d%d%2d: %6.2f -> %6.2f ",iSmType,iSm,iRpc,iBin,dCorP,dCorP+dYShift);
+ fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dYShift);
+ fhCalChannelDy->Fill(dYShift);
+ }
+ }
+ }
+ else { // no individuak y-corrections
+ for (Int_t iBin = 0; iBin < fhCorPos[iDetIndx]->GetNbinsX(); iBin++) {
+ Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
+ fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dYShift);
+ fhCalChannelDy->Fill(dYShift);
+ }
+ }
+ } break; // case iOpt0=0
+ case 1: // update channel mean
+ {
+ LOG(info) << "Update time offsets for TSR " << iSmType << iSm << iRpc << " with opt1 " << iOpt1;
TProfile* hpP = fhCalPos[iDetIndx]->ProfileX();
TProfile* hpT = fhCalTOff[iDetIndx]->ProfileX();
TH1* hCalP = fhCalPos[iDetIndx]->ProjectionX();
TH1* hCalT = fhCalTOff[iDetIndx]->ProjectionX();
//fhCorPos[iDetIndx]->Add((TH1 *)hpP,-1.);
//fhCorTOff[iDetIndx]->Add((TH1*)hpT,-1.);
- double dTOffMeanShift = 0.;
- double dNCh = 0.;
- for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
- Double_t dDt = hpT->GetBinContent(iBin + 1);
- Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
- Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
- Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
- Double_t dDp = hpP->GetBinContent(iBin + 1);
- Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
- LOG(info) << "Cts check for " << fhCalTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT << ", "
- << dCtsP << ", " << MINCTS;
- if (dCtsT > MINCTS && dCtsP > MINCTS) {
- // Fit Gaussian around peak
- TH1* hpPy =
- (TH1*) fhCalPos[iDetIndx]->ProjectionY(Form("PosPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
- Double_t dFMean = hpPy->GetBinCenter(hpPy->GetMaximumBin());
- Double_t dFLim = 0.5; // CAUTION, fixed numeric value
- Double_t dBinSize = hpPy->GetBinWidth(1);
- dFLim = TMath::Max(dFLim, 5. * dBinSize);
- TFitResultPtr fRes = hpPy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
- if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
- dDp = fRes->Parameter(1); //overwrite mean
- }
- TH1* hpTy =
- (TH1*) fhCalTOff[iDetIndx]->ProjectionY(Form("TOffPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
- dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
- dFLim = 0.5; // CAUTION, fixed numeric value
- dBinSize = hpTy->GetBinWidth(1);
+ if (iOpt1 == 3) { // update counter times from track pulls
+ TH1* hCalTpy = fhCalTOff[iDetIndx]->ProjectionY();
+ Double_t dCtsTpy = hCalTpy->GetEntries();
+ Double_t dDt = 0;
+ if (dCtsTpy > MINCTS) {
+ double dFMean = hCalTpy->GetBinCenter(hCalTpy->GetMaximumBin());
+ double dFLim = 0.5; // CAUTION, fixed numeric value
+ double dBinSize = hCalTpy->GetBinWidth(1);
dFLim = TMath::Max(dFLim, 5. * dBinSize);
- fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ TFitResultPtr fRes = hCalTpy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
- dDt = fRes->Parameter(1); //overwrite mean
- }
- // Double_t dDpRes = fRes->Parameter(2);
- if (iSmType == 5) { // do not shift beam counter in time
- fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt - dBeamTOff);
- dTOffMeanShift += dDt - dBeamTOff;
- }
- else {
- fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt + dBeamTOff);
- dTOffMeanShift += dDt + dBeamTOff;
- }
- dNCh += 1.;
- if (iOpt1 > 0) { // active y-position correction
- fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dDp);
- }
- if (iDetIndx > -1) {
- LOG(info) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, beam %6.3f, new %6.3f",
- fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dBeamTOff,
- dCorT + dDt + dBeamTOff);
+ dDt = fRes->Parameter(1);
+ LOG(info) << "Update cor hist " << fhCorTOff[iDetIndx]->GetName() << " by " << dDt << ", " << dBeamTOff;
+ if (iSmType == 5) { // do not shift beam counter in time
+ dDt -= dBeamTOff;
+ }
+ else {
+ dDt += dBeamTOff;
+ }
+ fhCalCounterDt->Fill(dDt);
+ for (int iCh = 0; iCh < fhCorTOff[iDetIndx]->GetNbinsX(); iCh++) {
+ fhCalChannelDt->Fill(dDt);
+ fhCorTOff[iDetIndx]->SetBinContent(iCh + 1, fhCorTOff[iDetIndx]->GetBinContent(iCh + 1) + dDt);
+ }
}
}
}
- if (dNCh > 0) dTOffMeanShift /= dNCh;
- LOG(info) << "Apply dTOffMeanShift " << dTOffMeanShift << " to " << fhCorTOff[iDetIndx]->GetName();
- for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) { //preserve mean offset
- fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, fhCorTOff[iDetIndx]->GetBinContent(iBin + 1) - dTOffMeanShift);
+ else {
+ double dTOffMeanShift = 0.;
+ double dNCh = 0.;
+ for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
+ Double_t dDt = hpT->GetBinContent(iBin + 1);
+ Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
+ Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
+ Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
+ Double_t dDp = hpP->GetBinContent(iBin + 1);
+ Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
+ LOG(debug) << "Cts check for " << fhCalTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
+ << ", " << dCtsP << ", " << MINCTS;
+ if (dCtsT > MINCTS && dCtsP > MINCTS) {
+ // Fit Gaussian around peak
+ TH1* hpPy =
+ (TH1*) fhCalPos[iDetIndx]->ProjectionY(Form("PosPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
+ if (hpPy->Integral() > MINCTS) {
+ Double_t dFMean = hpPy->GetBinCenter(hpPy->GetMaximumBin());
+ Double_t dFLim = 0.5; // CAUTION, fixed numeric value
+ Double_t dBinSize = hpPy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = hpPy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDp = fRes->Parameter(1); //overwrite mean
+ fhCalChannelDy->Fill(dDt);
+ }
+ }
+ TH1* hpTy =
+ (TH1*) fhCalTOff[iDetIndx]->ProjectionY(Form("TOffPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
+ if (hpTy->Integral() > MINCTS) {
+ Double_t dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
+ Double_t dFLim = 0.5; // CAUTION, fixed numeric value
+ Double_t dBinSize = hpTy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1); //overwrite mean
+ fhCalChannelDt->Fill(dDt);
+ }
+ }
+ // Double_t dDpRes = fRes->Parameter(2);
+ if (iSmType == 5) { // do not shift beam counter in time
+ fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt - dBeamTOff);
+ dTOffMeanShift += dDt - dBeamTOff;
+ }
+ else {
+ fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt + dBeamTOff);
+ dTOffMeanShift += dDt + dBeamTOff;
+ }
+ dNCh += 1.;
+ if (iOpt1 > 0) { // active y-position correction
+ fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dDp);
+ }
+ if (iDetIndx > -1) {
+ LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, beam %6.3f, new %6.3f",
+ fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dBeamTOff,
+ dCorT + dDt + dBeamTOff);
+ }
+ }
+ }
+ if (dNCh > 0) dTOffMeanShift /= dNCh;
+ //LOG(info) << "Apply dTOffMeanShift " << dTOffMeanShift << " to " << fhCorTOff[iDetIndx]->GetName();
+ for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) { //preserve mean offset
+ fhCorTOff[iDetIndx]->SetBinContent(iBin + 1,
+ fhCorTOff[iDetIndx]->GetBinContent(iBin + 1) - dTOffMeanShift);
+ }
}
} break;
- case 2: { // update individual channel walks
- LOG(info) << "Update Walks for TSR " << iSmType << iSm << iRpc;
- if (iSmType == 5) continue; // no walk correction for beam counter up to now
- const Double_t MinCounts = 10.;
- Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
- for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
- for (Int_t iSide = 0; iSide < 2; iSide++) {
- //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
- TProfile* hpW = fhCalWalk[iDetIndx][iCh][iSide]->ProfileX(); // mean deviation
- TH1* hCW = fhCalWalk[iDetIndx][iCh][iSide]->ProjectionX(); // contributing counts
- if (hCW->GetEntries() == 0) {
- LOG(info) << "No entries in " << hCW->GetName();
- continue;
- }
- Double_t dCorT = 0;
- for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
- Double_t dCts = hCW->GetBinContent(iBin + 1);
- Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1); // current value
- if (iBin > 0 && dCts == 0)
- fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1,
- fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin));
- if (dCts > MinCounts) { dCorT = hpW->GetBinContent(iBin + 1); }
- fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1, dWOff + dCorT); //set new value
+ case 2: { // iopt1==11: update individual channel walks from inside Cluster deviations
+ switch (iOpt1) {
+ case 10: {
+ LOG(debug) << "Update Cluster Offsets for TSR " << iSmType << iSm << iRpc;
+ if (iSmType == 5) continue; // do not shift beam counter in time
+ //TProfile* hpP = fhCalDelPos[iDetIndx]->ProfileX();
+ TProfile* hpT = fhCalDelTOff[iDetIndx]->ProfileX();
+ TH1* hCalP = fhCalDelPos[iDetIndx]->ProjectionX();
+ TH1* hCalT = fhCalDelTOff[iDetIndx]->ProjectionX();
+ //double dTOffMeanShift = 0.;
+ //double dNCh = 0.;
+ for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
+ Double_t dDt = hpT->GetBinContent(iBin + 1);
+ Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
+ Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
+ Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
+ //Double_t dDp = hpP->GetBinContent(iBin + 1);
+ //Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
+ LOG(debug) << "Cts check for " << fhCalDelTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
+ << ", " << dCtsP << ", " << MINCTS;
+ if (dCtsT > MINCTS) { //&& dCtsP > MINCTS) {
+ // Fit Gaussian around peak
+ TH1* hpTy = (TH1*) fhCalDelTOff[iDetIndx]->ProjectionY(Form("DelTOffPy_%d_%d", iDetIndx, iBin),
+ iBin + 1, iBin + 1);
+ double dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
+ double dFLim = 0.5; // CAUTION, fixed numeric value
+ double dBinSize = hpTy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1); //overwrite mean
+ }
+ fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt);
+ fhCalChannelDt->Fill(dDt);
+ }
+ if (iDetIndx > -1) {
+ LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, new %6.3f",
+ fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dCorT + dDt);
+ }
}
- // determine effective/count rate weighted mean
- Double_t dMean = 0;
- Double_t dCtsAll = 0.;
- for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
- Double_t dCts = hCW->GetBinContent(iBin + 1);
- Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1);
- if (dCts > MinCounts) {
- dCtsAll += dCts;
- dMean += dCts * dWOff;
+ } break;
+ case 12: {
+ LOG(debug) << "Update Cluster Position for TSR " << iSmType << iSm << iRpc;
+ if (iSmType == 5) continue; // do not shift beam counter in time
+ TProfile* hpP = fhCalDelPos[iDetIndx]->ProfileX();
+ TProfile* hpT = fhCalDelTOff[iDetIndx]->ProfileX();
+ TH1* hCalP = fhCalDelPos[iDetIndx]->ProjectionX();
+ TH1* hCalT = fhCalDelTOff[iDetIndx]->ProjectionX();
+ //double dTOffMeanShift = 0.;
+ //double dNCh = 0.;
+ for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
+ Double_t dDt = hpT->GetBinContent(iBin + 1);
+ Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
+ Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
+ Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
+ Double_t dDp = hpP->GetBinContent(iBin + 1);
+ Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
+ LOG(debug) << "Cts check for " << fhCalDelTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
+ << ", " << dCtsP << ", " << MINCTS;
+ if (dCtsT > MINCTS && dCtsP > MINCTS) {
+ // Fit Gaussian around peak
+ TH1* hpTy = (TH1*) fhCalDelTOff[iDetIndx]->ProjectionY(Form("DelTOffPy_%d_%d", iDetIndx, iBin),
+ iBin + 1, iBin + 1);
+ double dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
+ double dFLim = 0.5; // CAUTION, fixed numeric value
+ double dBinSize = hpTy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDt = fRes->Parameter(1); //overwrite mean
+ }
+ fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt);
+ fhCalChannelDt->Fill(dDt);
+
+ TH1* hpPy = (TH1*) fhCalDelPos[iDetIndx]->ProjectionY(Form("DelPosPy_%d_%d", iDetIndx, iBin),
+ iBin + 1, iBin + 1);
+ dFMean = hpPy->GetBinCenter(hpTy->GetMaximumBin());
+ dFLim = 0.5; // CAUTION, fixed numeric value
+ dBinSize = hpPy->GetBinWidth(1);
+ dFLim = TMath::Max(dFLim, 5. * dBinSize);
+ fRes = hpPy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ dDp = fRes->Parameter(1); //overwrite mean
+ }
+ fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dDp);
+ fhCalChannelDy->Fill(dDp);
+ }
+ if (iDetIndx > -1) {
+ LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, new %6.3f",
+ fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dCorT + dDt);
}
}
- if (dCtsAll > 0.) dMean /= dCtsAll;
-
- LOG(info) << "Mean shift for TSRCS " << iSmType << iSm << iRpc << iCh << iSide << ": " << dMean;
-
- // keep mean value at 0
- for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
- Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1); // current value
- fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1, dWOff - dMean); //set new value
- if (iCh == 0 && iSmType == 8 && iBin == 10) // debugging
- LOG(info) << "New Walk for " << fhCorWalk[iDetIndx][iCh][iSide]->GetName() << " bin " << iBin << ": "
- << fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1) << ", Mean " << dMean
- << ", WOff " << dWOff;
+ } break;
+ case 0:
+ case 11: {
+ LOG(debug) << "Update Walks for TSR " << iSmType << iSm << iRpc;
+ //if (iSmType == 5) continue; // no walk correction for beam counter up to now
+ const Double_t MinCounts = 10.;
+ Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
+ TH1* hCY = fhCalPos[iDetIndx]->ProjectionY(Form("%s_py_%d", fhCalPos[iDetIndx]->GetName(), iCh),
+ iCh + 1, iCh + 1);
+ if (hCY->GetEntries() > 1) fhCalChannelDy->Fill(hCY->GetRMS());
+ //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
+ TProfile* hpW0 = fhCalWalk[iDetIndx][iCh][0]->ProfileX(); // mean deviation
+ TH1* hCW0 = fhCalWalk[iDetIndx][iCh][0]->ProjectionX(); // contributing counts
+ TProfile* hpW1 = fhCalWalk[iDetIndx][iCh][1]->ProfileX(); // mean deviation
+ TH1* hCW1 = fhCalWalk[iDetIndx][iCh][1]->ProjectionX(); // contributing counts
+ if (iOpt0 == 3) { // use average distributions
+ hpW0 = fhCalWalkAv[iDetIndx]->ProfileX(); // mean deviation
+ hCW0 = fhCalWalkAv[iDetIndx]->ProjectionX(); // contributing counts
+ hpW1 = fhCalWalkAv[iDetIndx]->ProfileX(); // mean deviation
+ hCW1 = fhCalWalkAv[iDetIndx]->ProjectionX(); // contributing counts
+ }
+ if (hCW0->GetEntries() == 0 || hCW1->GetEntries() == 0) {
+ //LOG(info) << "No entries in " << hCW0->GetName();
+ continue;
+ }
+ Double_t dCorT = 0;
+ for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][0]->GetNbinsX(); iBin++) {
+ Double_t dCts0 = hCW0->GetBinContent(iBin + 1);
+ Double_t dCts1 = hCW1->GetBinContent(iBin + 1);
+ Double_t dWOff0 = fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin + 1); // current value
+ Double_t dWOff1 = fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin + 1); // current value
+ if (iBin > 0 && dCts0 == 0 && dCts1 == 0) {
+ fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1,
+ fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin));
+ fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1,
+ fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin));
+ }
+ dCorT = 0.;
+ if (dCts0 > MinCounts && dCts1 > MinCounts) {
+ dCorT = 0.5 * (hpW0->GetBinContent(iBin + 1) + hpW1->GetBinContent(iBin + 1));
+ fhCalChannelDt->Fill(dCorT);
+ }
+ fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1, dWOff0 + dCorT); //set new value
+ fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1, dWOff1 + dCorT); //set new value
+ if (iSmType == 0 && iSm == 0 && iRpc == 2) // debugging
+ LOG(info) << "UpdWalk " << fhCorWalk[iDetIndx][iCh][0]->GetName() << " bin " << iBin << " Tot "
+ << hCW0->GetBinCenter(iBin + 1) << ", " << hCW1->GetBinCenter(iBin + 1) << ": curVal "
+ << dWOff0 << ", " << dWOff1 << ", Cts " << dCts0 << ", " << dCts1 << ", dCorT " << dCorT
+ << ", newVal " << dWOff0 + dCorT << ", " << dWOff1 + dCorT;
+ }
+ // determine effective/count rate weighted mean
+ /*
+ Double_t dMean = 0;
+ Double_t dCtsAll = 0.;
+ for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
+ Double_t dCts = hCW->GetBinContent(iBin + 1);
+ Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1);
+ if (dCts > MinCounts) {
+ dCtsAll += dCts;
+ dMean += dCts * dWOff;
+ }
+ }
+ if (dCtsAll > 0.) dMean /= dCtsAll;
+
+ //LOG(info) << "Mean shift for TSRCS " << iSmType << iSm << iRpc << iCh << iSide << ": " << dMean;
+ // keep mean value at 0
+ for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
+ Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1); // current value
+ if (iSmType == 0 && iSm == 0 && iRpc == 2) // debugging
+ LOG(info) << "UpdWalk " << fhCorWalk[iDetIndx][iCh][iSide]->GetName()
+ << " bin " << iBin << " Tot " << hCW->GetBinCenter(iBin + 1)
+ << ": curDev " << hpW->GetBinContent(iBin + 1)
+ << ", Cts " << hCW->GetBinContent(iBin + 1)
+ << ", all " << dCtsAll
+ << ", oldWOff " << dWOff - hpW->GetBinContent(iBin + 1)
+ << ", Mean " << dMean
+ << ", newWOff " << dWOff - dMean;
+ fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1, dWOff + dMean); //set new value
+ }
+ */
}
- }
- }
+ } break; //iopt1=11 end
+ default:; LOG(info) << "Calibrator option " << iOpt << " not implemented ";
+ } // switch iOpt1 end
} break;
case 9: // update channel means on cluster level
@@ -772,7 +1828,8 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
if (NULL == fChannelInfo) LOG(fatal) << Form("invalid ChannelInfo for 0x%08x", iChId);
- double dThr = 10.;
+ double dThr = dCtsP / 100;
+ if (dThr < 3.) dThr = 3.;
double* dRes =
fTofClusterizer->find_yedges((const char*) (hPos_py->GetName()), dThr, fChannelInfo->GetSizey());
/*
@@ -791,9 +1848,11 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
}
}
- //TofOff
- TProfile* hpT = fhCalTofOff[iDetIndx]->ProfileX();
- TH1* hCalT = fhCalTofOff[iDetIndx]->ProjectionX();
+ //TofOff or TOff as input ?
+ //TProfile* hpT = fhCalTofOff[iDetIndx]->ProfileX();
+ //TH1* hCalT = fhCalTofOff[iDetIndx]->ProjectionX();
+ TProfile* hpT = fhCalTOff[iDetIndx]->ProfileX();
+ TH1* hCalT = fhCalTOff[iDetIndx]->ProjectionX();
for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
//Double_t dDt = hpT->GetBinContent(iBin + 1);
Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
@@ -812,16 +1871,20 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
//if (fRes == 0 )
if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
//if (TMath::Abs(TMean - fRes->Parameter(1)) > 1.)
+ /*
LOG(info) << "CalTOff "
<< Form("TSRC %d%d%d%d, %s gaus %8.2f %8.2f %8.2f for "
"TM %8.2f, TBeam %6.2f",
iSmType, iSm, iRpc, iBin, hTy->GetName(), fRes->Parameter(0), fRes->Parameter(1),
fRes->Parameter(2), dTmean, dBeamTOff);
+ */
dTmean = fRes->Parameter(1); //overwrite mean
}
}
+ /*
LOG(info) << Form("UpdateCalTOff TSRC %d%d%d%2d, cts %d: %6.2f -> %6.2f, %6.2f ", iSmType, iSm, iRpc,
iBin, (int) dCtsT, dCorT, dCorT + dTmean, dCorT + hpT->GetBinContent(iBin + 1));
+ */
fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dTmean);
}
}
@@ -834,31 +1897,41 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
double dTotMean = hbin->GetMean();
// Do gaus fit around maximum
Double_t dCtsTot = hbin->GetEntries();
- if (kFALSE && dCtsTot > MINCTS) {
+ if (dCtsTot > MINCTS) {
double dFMean = hbin->GetBinCenter(hbin->GetMaximumBin());
- double dFLim = dFMean; // CAUTION,
+ double dFLim = dFMean * 0.5; // CAUTION,
//LOG(info)<<Form("FitTot TSRC %d%d%d%2d: Mean %6.2f, Width %6.2f ",iSmType,iSm,iRpc,iBin,dFMean,dFLim);
if (dFMean > 2.) {
TFitResultPtr fRes = hbin->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
dTotMean = fRes->Parameter(1); //overwrite mean
}
+ else {
+ LOG(warn) << "FitFail " << hbin->GetName();
+ }
}
}
+ else { //append to broken channel list
+ int iCh = iBin % 2;
+ int iSide = iBin - iCh * 2;
+ CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, iCh, iSide);
+ Int_t iChId = fTofId->SetDetectorInfo(xDetInfo);
+ TString shcmd = Form("echo %d >> %s ", iChId, cBadChannelFile.Data());
+ gSystem->Exec(shcmd.Data());
+ }
if (dTotMean > 0.) { // prevent zero-divide
double dCorTot = fhCorTot[iDetIndx]->GetBinContent(ib);
/*
LOG(info)<<Form("UpdateCalTot TSRC %d%d%d%2d: %6.2f, %6.2f, %6.2f -> %6.2f ",iSmType,iSm,iRpc,iBin,
dCorTot,dCalTotMean,dTotMean,dCorTot*dTotMean/dCalTotMean);
- */
+ */
fhCorTot[iDetIndx]->SetBinContent(ib, dCorTot * dTotMean / dCalTotMean);
}
}
} break;
-
- default: LOG(warning) << "No valid calibration mode " << iOpt0;
- } //switch( iOpt) end
+ default: LOG(fatal) << "No valid calibration mode " << iOpt0;
+ } //switch( iOpt0) end
}
}
else {
@@ -870,6 +1943,9 @@ Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
WriteHist(fCalParFileNew);
fCalParFileNew->Close();
}
+ else {
+ WriteHist(fCalParFile);
+ }
fCalParFile->Close();
/// Restore old global file and folder pointer to avoid messing with FairRoot
@@ -921,11 +1997,19 @@ void CbmTofCalibrator::ReadHist(TFile* fHist)
fhCorWalk[iDetIndx][iCh].resize(2);
for (Int_t iSide = 0; iSide < 2; iSide++) {
//LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
- fhCorWalk[iDetIndx][iCh][iSide] = (TH1*) gDirectory->FindObjectAny(
- Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_Walk_px", iSmType, iSm, iRpc, iCh, iSide));
+ TString hname = Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_Walk_px", iSmType, iSm, iRpc, iCh, iSide);
+ fhCorWalk[iDetIndx][iCh][iSide] = (TH1*) gDirectory->FindObjectAny(hname);
if (NULL == fhCorWalk[iDetIndx][iCh][iSide]) {
LOG(warn) << "No Walk histo for TSRCS " << iSmType << iSm << iRpc << iCh << iSide;
- continue;
+ if (NULL != fhCalWalk[iDetIndx][iCh][iSide]) {
+ LOG(info) << "Create walk correction histo " << hname;
+ fhCorWalk[iDetIndx][iCh][iSide] =
+ new TH1D(hname,
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSm,
+ iRpc, iCh, iSide),
+ nbClWalkBinX, fTofClusterizer->GetTotMin(), fTofClusterizer->GetTotMax());
+ }
+ //continue;
}
if (iSmType == 8 && iSide == 1) { //pad special treatment
@@ -938,6 +2022,35 @@ void CbmTofCalibrator::ReadHist(TFile* fHist)
}
}
}
+
+ int iNbTypes = fDigiBdfPar->GetNbSmTypes();
+ hSvel.resize(iNbTypes);
+ for (int iSmType = 0; iSmType < iNbTypes; iSmType++) {
+ int iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ if (iNbRpc > 0 && fDigiBdfPar->GetNbSm(iSmType) > 0) {
+ TH1* hTmp = (TH1*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
+ if (hTmp == NULL) {
+ TDirectory* oldir = gDirectory;
+ gROOT->cd();
+ LOG(warn) << Form("cl_SmT%01d_Svel not found, creating ... in %s ", iSmType, gDirectory->GetName());
+ hSvel[iSmType] =
+ new TH1F(Form("cl_SmT%01d_Svel", iSmType), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iSmType),
+ fDigiBdfPar->GetNbSm(iSmType) * iNbRpc, 0, fDigiBdfPar->GetNbSm(iSmType) * iNbRpc);
+
+ for (int i = 0; i < hSvel[iSmType]->GetNbinsX(); i++)
+ hSvel[iSmType]->SetBinContent(i + 1, 1.); // set default
+ oldir->cd();
+ }
+ else {
+ //LOG(info) << "hSvel "<< hTmp->GetName() << " located in " << gDirectory->GetName();
+ TDirectory* oldir = gDirectory;
+ gROOT->cd();
+ hSvel[iSmType] = (TProfile*) hTmp->Clone();
+ //LOG(info) << "hSvel "<< hSvel[iSmType]->GetName() << " cloned in " << gDirectory->GetName();
+ oldir->cd();
+ }
+ }
+ }
}
void CbmTofCalibrator::WriteHist(TFile* fHist)
@@ -955,10 +2068,23 @@ void CbmTofCalibrator::WriteHist(TFile* fHist)
Int_t iNbCh = (Int_t) fhCorWalk[iDetIndx].size();
for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
for (Int_t iSide = 0; iSide < 2; iSide++) {
- fhCorWalk[iDetIndx][iCh][iSide]->Write();
+ //LOG(info)<<"Write " << fhCorWalk[iDetIndx][iCh][iSide]->GetName();
+ if (NULL != fhCorWalk[iDetIndx][iCh][iSide]) fhCorWalk[iDetIndx][iCh][iSide]->Write();
+ else
+ LOG(warn) << "Walk correction histo for DetIndx " << iDetIndx << ", Ch " << iCh << ", S " << iSide
+ << " not found";
}
}
}
+ for (Int_t iS = 0; iS < fDigiBdfPar->GetNbSmTypes(); iS++) {
+ if (NULL != hSvel[iS]) {
+ hSvel[iS]->Write();
+ LOG(warn) << "Wrote " << hSvel[iS]->GetName() << " to " << gDirectory->GetName();
+ }
+ else {
+ LOG(warn) << Form("cl_SmT%01d_Svel not found ", iS);
+ }
+ }
oldir->cd();
}
@@ -1013,5 +2139,269 @@ void CbmTofCalibrator::HstDoublets(CbmTofTracklet* pTrk)
}
}
+double* CbmTofCalibrator::find_tofedge(const char* hname, double dThr, double dLen)
+{
+ TH1* h1 = (TH1*) gROOT->FindObjectAny(hname);
+ static double dRes[2] = {2 * 0};
+ dRes[0] = 0.;
+ if (NULL != h1) {
+ /*
+ LOG(debug) << "Inspect " << h1->GetName() << ", entries " << h1->GetNbinsX()
+ << ", Thr " << dThr << ", Len " << dLen << endl;
+ */
+
+ const int iMaxInt = 3;
+ int iMax = h1->GetMaximumBin();
+ int iLow = TMath::Max(iMax - iMaxInt, 1);
+ int iHigh = TMath::Min(iMax + iMaxInt, h1->GetNbinsX());
+ double dMax = 0.;
+ int nMax = 0;
+ for (int i = iLow; i < iHigh; i++) {
+ dMax += h1->GetBinContent(i);
+ nMax++;
+ }
+ dMax /= nMax;
+
+ // determine average noise and RMS per bin
+ double dNoise = 0;
+ double dNoise2 = 0;
+ double dRMS = 0;
+ int iNbins = 0;
+ for (int iBin = 0; iBin < h1->GetNbinsX() / 4; iBin++) {
+ if (h1->GetBinContent(iBin + 1) > 0 || iNbins > 0) { // start counting from first bin with entries
+ dNoise += h1->GetBinContent(iBin + 1);
+ dNoise2 += h1->GetBinContent(iBin + 1) * h1->GetBinContent(iBin + 1);
+ iNbins++;
+ }
+ }
+ if (iNbins > 0) {
+ dNoise /= iNbins;
+ dNoise2 /= iNbins;
+ dRMS = TMath::Sqrt(dNoise2 - dNoise * dNoise);
+ }
+ else {
+ dNoise = h1->GetMaximum() * 0.1;
+ LOG(warn) << h1->GetName() << ": Noise level could not be determined, init to " << dNoise;
+ }
+ double dLev = dNoise + (dMax - dNoise) * dThr;
+ dLev = TMath::Max(dNoise + 3. * dRMS, dLev);
+ if (dLev == 0) {
+ LOG(warn) << "Invalid threshold level " << dLev;
+ return dRes;
+ }
+ int iBl = -1;
+ double xLow = h1->GetBinCenter(1);
+
+ for (int iBin = 1; iBin < h1->GetNbinsX() - dLen; iBin++) {
+ int iLen = 0;
+ for (; iLen < dLen; iLen++) {
+ LOG(debug) << h1->GetName() << ": Lev " << dLev << ", Bin " << iBin << ", Len " << iLen << ", cont "
+ << h1->GetBinContent(iBin + iLen);
+ if (h1->GetBinContent(iBin + iLen) < dLev) break;
+ }
+ if (iLen == dLen) {
+ if (iBin == 1) {
+ LOG(warn) << "find_tofedge: " << h1->GetName() << ", Lvl1 " << dLev;
+ iBin--;
+ dLev *= 2;
+ continue;
+ }
+ iBl = iBin;
+ xLow = h1->GetBinCenter(iBin);
+ if (iBin == 1) { // not active since did not work, case already caught above
+ LOG(warn) << GetName() << ": " << h1->GetName() << " Lvl1 " << dLev;
+ xLow *= 2;
+ }
+ break;
+ }
+ }
+ dRes[0] = xLow;
+ LOG(info) << h1->GetName() << " with Thr " << dLev << ", Noise " << dNoise << ", RMS " << dRMS << " at iBl " << iBl
+ << ", TOff " << dRes[0];
+ }
+ return dRes;
+}
+
+double* CbmTofCalibrator::find_tofedge(const char* hname)
+{
+ double dThr = 0.1;
+ double dLen = 3.;
+ return find_tofedge(hname, dThr, dLen);
+}
+
+double CbmTofCalibrator::CalcChi2(TH1* h1, TH1* h2, int i)
+{
+ double chi2 = 0.;
+ double nBin = 0.;
+ double nEntries1 = h1->GetEntries();
+ double nEntries2 = h2->GetEntries();
+ const int CompRange = (int) h1->GetNbinsX() / 4;
+ for (int iBin1 = CompRange; iBin1 < h1->GetNbinsX() - CompRange; iBin1++) {
+ double diff = h1->GetBinContent(iBin1) / nEntries1 - h2->GetBinContent(iBin1 + i) / nEntries2;
+ chi2 += diff * diff;
+ nBin++;
+ }
+ chi2 /= nBin;
+ return chi2;
+}
+
+
+Double_t CbmTofCalibrator::f1_tedge(double* x, double* par)
+{
+ double xx = x[0];
+ //double wx = par[0] + par[1] * TMath::Exp(xx*par[2]);
+ double wx = par[0] + par[1] * (1. + TMath::Erf((xx - par[2]) / par[3]));
+
+ return wx;
+}
+
+double* CbmTofCalibrator::fit_tofedge(const char* hname, Double_t TMax, Double_t dThr)
+{
+ TH1* h1;
+ static double res[10] = {10 * 0.};
+ double err[10];
+ h1 = (TH1*) gROOT->FindObjectAny(hname);
+ if (NULL != h1) {
+ const double MinCts = 1000.;
+ res[2] = 0.;
+ if (h1->GetEntries() < MinCts) {
+ LOG(warn) << h1->GetName() << ": too few entries for edgefit " << h1->GetEntries();
+ return &res[2];
+ }
+
+ TAxis* xaxis = h1->GetXaxis();
+ Double_t Tmin = xaxis->GetXmin();
+ Double_t Tmax = xaxis->GetXmax();
+ Double_t dXmax = Tmax;
+ if (TMax == 0.) TMax = Tmax;
+ else
+ Tmax = TMax - 1.;
+ Tmax = TMath::Max(Tmax, Tmin + 0.3);
+ TFitResultPtr fRes = h1->Fit("pol0", "SQM", "", Tmin, Tmax);
+ TF1* f0 = h1->GetFunction("pol0");
+ if (fRes != 0 && f0 != NULL && (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED"))) {
+ //double dBackground=fRes->Parameter(0) ;
+ double dBackground = f0->GetParameter(0);
+ //cout << h1->GetName() << ": av background rate \t" << dBackground<< endl;
+ TF1* f1 = new TF1("TEdge", f1_tedge, Tmin, dXmax, 4);
+ if (TMax != 0.) {
+ f1->SetParameters(dBackground, 1000., 0., 0.1);
+ f1->SetParLimits(0, dBackground, dBackground); //fix background
+ }
+ f1->SetParLimits(2, -0.5, 0.5);
+ f1->SetParLimits(3, 0.05, 0.5);
+
+ h1->Fit("TEdge", "SQM", "", Tmin, TMax);
+ res[9] = f1->GetChisquare();
+
+ for (int i = 0; i < 4; i++) {
+ res[i] = f1->GetParameter(i);
+ err[i] = f1->GetParError(i);
+ }
+
+ if (res[1] < 0.) {
+ LOG(warn) << h1->GetName() << "fitted wrong sign of step function ";
+ res[0] = h1->GetBinCenter(h1->GetNbinsX());
+ return &res[0];
+ }
+ // determine toff at 10*background
+ /*
+ double dToff=0.;
+ if (res[1] !=0. && res[2] != 0.) dToff=TMath::Log(10*dBackground/res[1])/res[2];
+ LOG(info) << "TEdge Fit of " << hname
+ << " in [" << Tmin <<","<< TMax <<"]"
+ << " ended with chi2 = " << res[9]
+ << ", Toff " << dToff
+ << Form(", noise level %7.2f +/- %5.2f, e0 %7.2f "
+ "+/- %5.2f, exp = %7.2f +/- %5.2f",
+ res[0], err[0], res[1], err[1], res[2], err[2]);
+ res[0]=dToff;
+ */
+ // FIXME: code duplication from find_tofedge
+ const int iMaxInt = 3;
+ int iMax = h1->GetMaximumBin();
+ int iLow = TMath::Max(iMax - iMaxInt, 1);
+ int iHigh = TMath::Min(iMax + iMaxInt, h1->GetNbinsX());
+ double dMax = 0.;
+ int nMax = 0;
+ for (int i = iLow; i < iHigh; i++) {
+ dMax += h1->GetBinContent(i);
+ nMax++;
+ }
+ dMax /= nMax;
+ double dLev = dThr * dMax;
+ dLev = dBackground + (dMax - dBackground) * 0.1;
+
+ double dToff = f1->GetX(dLev, Tmin, dXmax);
+
+ LOG(info) << "TEdge Fit of " << hname << " in [" << Tmin << "," << TMax << "/" << dXmax << "]"
+ << ": chi2 " << res[9]
+ << Form(
+ ", noise %6.2f +/-%4.2f, norm %6.2f +/-%4.2f, mean %6.2f +/-%5.3f, wid %6.2f +/-%5.3f -> %6.3f ",
+ res[0], err[0], res[1], err[1], res[2], err[2], res[3], err[3], dToff);
+ res[2] = dToff;
+ }
+ else {
+ LOG(warn) << "pol0 fit failed, Minuit " << gMinuit->fCstatu << ", fRes " << fRes << ", "
+ << ", f0 " << f0;
+ }
+ }
+ return &res[2];
+}
+
+double* CbmTofCalibrator::fit_tofedge(const char* hname)
+{
+ Double_t Tmax = 0.3;
+ return fit_tofedge(hname, Tmax, fTofClusterizer->GetEdgeThr());
+}
+
+double CbmTofCalibrator::TruncatedMeanY(TH2* h2, double dRmsLim)
+{
+ double dMeanIn = h2->ProjectionY()->GetMean();
+ double dRmsIn = h2->ProjectionY()->GetRMS();
+ const int Nx = h2->GetNbinsX();
+
+ int Ntru = 0;
+ TH1* h1[Nx];
+ double dMean = 0.;
+ double dMean2 = 0.;
+ double dNorm = 0.;
+ double dRms = 0.;
+ for (int i = 0; i < Nx; i++) {
+ h1[i] = h2->ProjectionY(Form("%s_py%d", h2->GetName(), i), i + 1, i + 1);
+ if (TMath::Abs(dMeanIn - h1[i]->GetMean()) < dRmsLim * dRmsIn) {
+ dMean += h1[i]->GetMean() * h1[i]->GetEntries();
+ dMean2 += h1[i]->GetMean() * h1[i]->GetMean() * h1[i]->GetEntries();
+ dNorm += h1[i]->GetEntries();
+ Ntru++;
+ }
+ }
+ if (Ntru > 0) {
+ dMean /= dNorm;
+ dMean2 /= dNorm;
+ if (dMean2 > dMean * dMean) { dRms = TMath::Sqrt(dMean2 - dMean * dMean); }
+ else {
+ LOG(error) << h2->GetName() << " Invalid Rms calculation: " << dMean2 << ", " << dMean * dMean;
+ return dMean;
+ }
+ double dMean1 = dMean;
+ dMean = 0.;
+ dMean2 = 0.;
+ dNorm = 0.;
+ Ntru = 0;
+ for (int j = 0; j < Nx; j++) {
+ if (TMath::Abs(dMean1 - h1[j]->GetMean()) < dRmsLim * dRms) {
+ dMean += h1[j]->GetMean() * h1[j]->GetEntries();
+ dMean2 += h1[j]->GetMean() * h1[j]->GetMean() * h1[j]->GetEntries();
+ dNorm += h1[j]->GetEntries();
+ Ntru++;
+ }
+ }
+ if (Ntru > 0) { dMean /= dNorm; }
+ }
+ LOG(info) << h2->GetName() << ": TruncY " << dMean << "( " << Ntru << ")"
+ << ", full " << dMeanIn << ", RMS " << dRms;
+ return dMean;
+}
ClassImp(CbmTofCalibrator)
diff --git a/reco/detectors/tof/CbmTofCalibrator.h b/reco/detectors/tof/CbmTofCalibrator.h
index f89430cbac6e5f1ead5523e88630b8861b4ce409..c90c5aa5f1ce382d59ac9fc43c119c2d6fde1d7a 100644
--- a/reco/detectors/tof/CbmTofCalibrator.h
+++ b/reco/detectors/tof/CbmTofCalibrator.h
@@ -63,11 +63,19 @@ public:
Bool_t InitParameters();
Bool_t CreateCalHist();
void FillCalHist(CbmTofHit* pHit, Int_t iOpt, CbmEvent* pEvent = NULL);
+ void FillHitCalHist(CbmTofHit* pHit, Int_t iOpt, CbmEvent* pEvent = NULL, TClonesArray* tHitColl = NULL);
void FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* pEvent = NULL);
Bool_t UpdateCalHist(Int_t iOpt);
void ReadHist(TFile* fhFile);
void WriteHist(TFile* fhFile);
void HstDoublets(CbmTofTracklet* pTrk);
+ double* find_tofedge(const char* hname, Double_t dThr, Double_t dLen);
+ double* find_tofedge(const char* hname);
+ double CalcChi2(TH1* h1, TH1* h2, int iShift);
+ double* fit_tofedge(const char* hname, Double_t dTmax, Double_t dThr);
+ double* fit_tofedge(const char* hname);
+ static double f1_tedge(double* x, double* par);
+ double TruncatedMeanY(TH2* pHst, double RmsLim = 1.);
inline void SetR0Lim(Double_t dVal) { fdR0Lim = dVal; }
inline void SetBeam(Bool_t bVal) { fbBeam = bVal; }
@@ -82,6 +90,8 @@ private:
CbmTofDigiPar* fDigiPar;
CbmTofDigiBdfPar* fDigiBdfPar;
+ CbmTofGeoHandler* fGeoHandler;
+ CbmTofDetectorId* fTofId;
TClonesArray* fTofDigiMatchColl; // TOF Digi Links
FairEventHeader* fEvtHeader;
@@ -90,12 +100,27 @@ private:
TH1* fhCalDX0;
TH1* fhCalDY0;
+ TH1* fhCalCounterDt;
+ TH1* fhCalCounterDy;
+ TH1* fhCalChannelDt;
+ TH1* fhCalChannelDy;
+
std::vector<TH2*> fhCalTot; // [nbDet]
std::vector<TH2*> fhCalPosition; // [nbDet]
std::vector<TH2*> fhCalPos; // [nbDet]
std::vector<TH2*> fhCalTOff; // [nbDet]
std::vector<TH2*> fhCalTofOff; // [nbDet]
+ std::vector<TH2*> fhCalDelPos; // [nbDet]
+ std::vector<TH2*> fhCalDelTOff; // [nbDet]
+ std::vector<TH2*> fhCalCluTrms; // [nbDet]
+ std::vector<TH2*> fhCalCluSize; // [nbDet]
+ std::vector<TH2*> fhCalWalkAv; // [nbDet]
std::vector<std::vector<std::vector<TH2*>>> fhCalWalk; // [nbDet][nbCh][nSide]
+ std::vector<std::vector<std::vector<TH2*>>> fhCalDtWalk; // [nbDet][nbCh][nSide]
+ std::vector<TH3*> fhCalXYTOff; // [nbDet]
+ std::vector<TH3*> fhCalXYTot; // [nbDet]
+ std::vector<std::vector<std::vector<TH3*>>> fhCalTotYWalk; // [nbDet][nbCh][nSide]
+ std::vector<std::vector<std::vector<TH3*>>> fhCalTotYTOff; // [nbDet][nbCh][nSide]
std::vector<TH1*> fhCorPos; // [nbDet]
std::vector<TH1*> fhCorTOff; // [nbDet]
diff --git a/reco/detectors/tof/CbmTofEventClusterizer.cxx b/reco/detectors/tof/CbmTofEventClusterizer.cxx
index 9d9398357cfd3942e802e025368d5f23096e5d26..911f6ccf8c791ee820dfdf5ff164e533d888415c 100644
--- a/reco/detectors/tof/CbmTofEventClusterizer.cxx
+++ b/reco/detectors/tof/CbmTofEventClusterizer.cxx
@@ -19,7 +19,8 @@
#include "CbmDigiManager.h"
#include "CbmEvent.h"
#include "CbmMatch.h"
-#include "CbmTofAddress.h" // in cbmdata/tof
+#include "CbmTimeSlice.h"
+#include "CbmTofAddress.h" // in cbmdata/tof
#include "CbmTofCalibrator.h"
#include "CbmTofCell.h" // in tof/TofData
#include "CbmTofCreateDigiPar.h" // in tof/TofTools
@@ -32,8 +33,11 @@
#include "CbmTofGeoHandler.h" // in tof/TofTools
#include "CbmTofHit.h" // in cbmdata/tof
#include "CbmTofPoint.h" // in cbmdata/tof
+#include "CbmTsEventHeader.h"
#include "CbmVertex.h"
+#include "TimesliceMetaData.h"
+
#include "TTrbHeader.h"
// CBMroot classes and includes
@@ -55,6 +59,7 @@
#include "TFitResult.h"
#include "TFitter.h"
#include "TGeoManager.h"
+#include "TGeoPhysicalNode.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
@@ -85,13 +90,14 @@ static Int_t iNSpill = 0;
static Int_t iNbTs = 0;
static int fiHitStart = 0;
static double dTsStartTime = 0.;
+static double dTsStartTimeLast = -1.;
const Double_t fdSpillDuration = 2.; // in seconds
const Double_t fdSpillBreak = 0.9; // in seconds
const Double_t dTimeRes = 0.08; // in ns
static Bool_t bAddBeamCounterSideDigi = kTRUE;
static TRandom3* fRndm = new TRandom3();
-
+std::vector<CbmTofDigi>* fT0DigiVec; //! T0 Digis
// std::vector< CbmTofPoint* > vPtsRef;
CbmTofEventClusterizer* CbmTofEventClusterizer::fInstance = 0;
@@ -111,6 +117,8 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fDigiBdfPar(NULL)
, fTrbHeader(NULL)
, fEvtHeader(NULL)
+ , fTsHeader(NULL)
+ , fTimeSlice(NULL)
, fTofDigiPointMatches(NULL)
, fDigiMan(nullptr)
, fEventsColl(nullptr)
@@ -136,6 +144,7 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fvdDifCh()
, fTofCalibrator(NULL)
, fhClustBuildTime(NULL)
+ , fhClustHitsDigi(NULL)
, fhHitsPerTracks(NULL)
, fhPtsPerHit(NULL)
, fhTimeResSingHits(NULL)
@@ -165,6 +174,10 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fhRpcDigiDTFD()
, fhRpcDigiDTMul()
, fhRpcDigiRate()
+ , fhRpcDigiTotLeft()
+ , fhRpcDigiTotRight()
+ , fhRpcDigiTotDiff()
+ , fhRpcDigiTotMap()
, fhRpcCluMul()
, fhRpcCluRate()
, fhRpcCluRate10s()
@@ -216,6 +229,9 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fvCPTotGain()
, fvCPTotOff()
, fvCPWalk()
+ , fvCPTOffY()
+ , fvCPTOffYBinWidth()
+ , fvCPTOffYRange()
, fvLastHits()
, fvDeadStrips()
, fvTimeLastDigi()
@@ -253,7 +269,10 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fSel2Rpc(0)
, fSel2Addr(0)
, fSel2MulMax(1)
+ , fiCluSizeMin(0)
+ , fNbCalHitMin(0)
, fDetIdIndexMap()
+ , fCellIdGeoMap()
, fviDetId()
, fPosYMaxScal(1.)
, fTRefDifMax(10.)
@@ -281,6 +300,11 @@ CbmTofEventClusterizer::CbmTofEventClusterizer(const char* name, Int_t verbose,
, fdTTotMean(2.)
, fdMaxTimeDist(0.)
, fdMaxSpaceDist(0.)
+ , fdEdgeThr(0.1)
+ , fdEdgeLen(3.)
+ , fdEdgeTbias(-0.3)
+ , fdEdgeFrange(0.4)
+ , fdModifySigvel(1.0)
, fdEvent(0)
, fdStartAnalysisTime(0.)
, fbSwapChannelSides(kFALSE)
@@ -316,7 +340,10 @@ InitStatus CbmTofEventClusterizer::Init()
if (kFALSE == CreateHistos()) return kFATAL;
- if (fCalMode % 10 > 7 && fDutId > -1) {
+ if ((fCalMode % 10 > 7 && fDutId > -1) || (fCalMode > 99 && fDutId < 0)) {
+ if (fiBeamRefType > -1)
+ fiBeamRefAddr = CbmTofAddress::GetUniqueAddress(fiBeamRefSm, fiBeamRefDet, 0, 0, fiBeamRefType);
+ LOG(info) << "Initialize Calibrator for Clusterizer ";
fTofCalibrator = new CbmTofCalibrator();
if (fTofCalibrator->Init() != kSUCCESS) return kFATAL;
return kSUCCESS;
@@ -373,9 +400,10 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
if (fTofCalDigiVecOut) fTofCalDigiVecOut->clear();
if (fEventsColl) {
- LOG(info) << "New timeslice " << iNbTs << " with " << fEventsColl->GetEntriesFast() << " events, "
- << fDigiMan->GetNofDigis(ECbmModuleId::kTof) << " TOF digis + "
- << fDigiMan->GetNofDigis(ECbmModuleId::kT0) << " T0 digis ";
+ if (NULL != fTsHeader)
+ LOG(info) << "New Ts " << iNbTs << ", size " << fEventsColl->GetSize() << " at " << fTsHeader->GetTsStartTime()
+ << " with " << fEventsColl->GetEntriesFast() << " events, " << fDigiMan->GetNofDigis(ECbmModuleId::kTof)
+ << " TOF digis + " << fDigiMan->GetNofDigis(ECbmModuleId::kT0) << " T0 digis ";
TStopwatch timerTs;
timerTs.Start();
@@ -387,6 +415,7 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
Int_t iNbHits = 0;
Int_t iNbCalDigis = 0;
fTofDigiMatchCollOut->Delete(); // costly, FIXME
+ //for(int i=0; i<fTofHitsCollOut->GetEntriesFast(); i++) delete fTofHitsCollOut->At(i);
fTofHitsCollOut->Delete(); // costly, FIXME
//fTofDigiMatchCollOut->Clear("C"); // not sufficient, memory leak
for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
@@ -400,8 +429,10 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
for (Int_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kT0Digi); iDigi++) {
Int_t iDigiIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kT0Digi, iDigi));
CbmTofDigi tDigi(fDigiMan->Get<CbmBmonDigi>(iDigiIndex));
- if (tDigi.GetType() != 5)
- LOG(fatal) << "Wrong T0 type " << tDigi.GetType() << ", Addr 0x" << std::hex << tDigi.GetAddress();
+ if (tDigi.GetType() != 5) {
+ //LOG(fatal) << "Wrong T0 type " << tDigi.GetType() << ", Addr 0x" << std::hex << tDigi.GetAddress();
+ tDigi.SetAddress(0, 0, 0, 0, 5); // convert to Tof Address
+ }
if (tDigi.GetSide() == 1) { // HACK for May22 setup
tDigi.SetAddress(tDigi.GetSm(), tDigi.GetRpc(), tDigi.GetChannel() + 6, 0, tDigi.GetType());
}
@@ -459,12 +490,19 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
//tEvent->SetMatch( (CbmMatch*)((*fTofDigiMatchCollOut)[iNbHits]) );
iNbHits++;
}
- if (fCalMode % 10 == 9 && fDutId > -1) {
- if (iNbHits - iHit0 > 0) { // outsource all calibration actions
+
+ if ((fCalMode % 10 == 9 && fDutId > -1) || (fCalMode > 99 && fDutId < 0)) {
+ if (iNbHits - iHit0 > fNbCalHitMin) { // outsource all calibration actions
LOG(debug) << "Call Calibrator with fiHitStart = " << fiHitStart << ", entries "
<< fTofHitsCollOut->GetEntriesFast();
CbmTofHit* pHit = (CbmTofHit*) fTofHitsCollOut->At(fiHitStart); // use most early hit as reference
- fTofCalibrator->FillCalHist(pHit, fCalMode, tEvent);
+ if (fCalMode < 100) fTofCalibrator->FillCalHist(pHit, fCalMode, tEvent);
+ else {
+ if (fCalMode < 1000) {
+ CbmTofHit* pRef = NULL;
+ fTofCalibrator->FillHitCalHist(pRef, fCalMode, tEvent, fTofHitsCollOut);
+ }
+ }
}
}
@@ -473,7 +511,7 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
fTofDigiVec.clear();
//fTofCalDigi->Delete();//Clear("C"); //otherwise memoryleak! FIXME
fTofCalDigiVec->clear();
- fTofHitsColl->Clear("C");
+ fTofHitsColl->Delete(); //Clear("C");
fTofDigiMatchColl->Delete(); //Clear("C");
}
if (0 < fiNbSkip1) {
@@ -512,6 +550,23 @@ void CbmTofEventClusterizer::Exec(Option_t* option)
fTofDigiVec.clear();
//if (fTofDigisColl) fTofDigisColl->Clear("C");
// Int_t iNbDigis=0; (VF) not used
+ if (NULL != fT0DigiVec) // 2022 data
+ for (Int_t iDigi = 0; iDigi < (int) (fT0DigiVec->size()); iDigi++) {
+ CbmTofDigi tDigi = fT0DigiVec->at(iDigi);
+ if (tDigi.GetType() != 5)
+ LOG(fatal) << "Wrong T0 type " << tDigi.GetType() << ", Addr 0x" << std::hex << tDigi.GetAddress();
+ if (tDigi.GetSide() == 1) { // HACK for May22 setup
+ tDigi.SetAddress(tDigi.GetSm(), tDigi.GetRpc(), tDigi.GetChannel() + 6, 0, tDigi.GetType());
+ }
+ fTofDigiVec.push_back(CbmTofDigi(tDigi));
+ }
+ else { // MC
+ /*
+ for (Int_t iDigi = 0; iDigi < fDigiMan->GetNofDigis(ECbmModuleId::kT0); iDigi++) {
+ const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigi);
+ }
+ */
+ }
for (Int_t iDigi = 0; iDigi < fDigiMan->GetNofDigis(ECbmModuleId::kTof); iDigi++) {
const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigi);
fTofDigiVec.push_back(CbmTofDigi(*tDigi));
@@ -526,8 +581,8 @@ void CbmTofEventClusterizer::ExecEvent(Option_t* /*option*/, CbmEvent* tEvent)
// Clear output arrays
//fTofCalDigisColl->Delete(); //otherwise memoryleak if 'CbmDigi::fMatch' points to valid MC match objects (simulation)! FIXME
fTofCalDigiVec->clear();
- fTofHitsColl->Clear("C");
- //fTofHitsColl->Delete(); // Computationally costly!, but hopefully safe
+ //fTofHitsColl->Clear("C");
+ fTofHitsColl->Delete(); // Computationally costly!, but hopefully safe
//for (Int_t i=0; i<fTofDigiMatchColl->GetEntriesFast(); i++) ((CbmMatch *)(fTofDigiMatchColl->At(i)))->ClearLinks(); // FIXME, try to tamper memory leak (did not help)
//fTofDigiMatchColl->Clear("C+L"); // leads to memory leak
fTofDigiMatchColl->Delete();
@@ -535,11 +590,21 @@ void CbmTofEventClusterizer::ExecEvent(Option_t* /*option*/, CbmEvent* tEvent)
if (bla) fiOutputTreeEntry = ((FairRootFileSink*) FairRootManager::Instance()->GetSink())->GetOutTree()->GetEntries();
// Check for corruption
- if (fTofDigiVec.size() > 10. * fDigiBdfPar->GetNbDet()) return; // FIXME constant in code, skip this event
-
+ if (fTofDigiVec.size() > 20. * fDigiBdfPar->GetNbDet()) { // FIXME constant in code, skip this event
+ LOG(warn) << "Skip processing of Tof DigiEvent with " << fTofDigiVec.size() << " digis ";
+ return;
+ }
fiNbHits = 0;
- if (fEvtHeader != NULL) dTsStartTime = fEvtHeader->GetEventTime();
- LOG(debug) << Form("BuildClu for event %d at tTS = %f ", (int) fdEvent, dTsStartTime);
+ if (fEvtHeader != NULL) {
+ dTsStartTime = fEvtHeader->GetEventTime();
+ if (fTsHeader != NULL) dTsStartTime = static_cast<double>(fTsHeader->GetTsStartTime());
+ if (dTsStartTime != dTsStartTimeLast) {
+ LOG(debug) << Form("BuildClu for event %d at tTS = %f ", (int) fdEvent, dTsStartTime);
+ dTsStartTimeLast = dTsStartTime;
+ }
+ }
+ else
+ LOG(fatal) << "No EvtHeader found!";
fStart.Set();
@@ -551,13 +616,12 @@ void CbmTofEventClusterizer::ExecEvent(Option_t* /*option*/, CbmEvent* tEvent)
fdEvent++;
FillHistos(tEvent);
-
- // fTofDigisColl->RemoveAll();
}
/************************************************************************************/
void CbmTofEventClusterizer::Finish()
{
+ FairLogger::GetLogger()->SetLogScreenLevel("info");
LOG(info) << "Finish with " << fdEvent << " processed events";
/// PAL: add run statistics for monitoring and perf evaluation
@@ -575,12 +639,14 @@ void CbmTofEventClusterizer::Finish()
LOG(info) << "=====================================";
if (fdEvent < 100) return; // don't save histos with insufficient statistics
- if (fDutId < 0) return; // no histograms were produced
- if (fCalMode % 10 < 8) { WriteHistos(); }
+ if (fCalMode % 10 < 8 && fDutId > -1) { WriteHistos(); }
else {
- fTofCalibrator->UpdateCalHist(fCalMode);
+ if ((fCalMode % 10 > 7 && fDutId > -1) || (fCalMode > 99 && fDutId < 0)) {
+ fTofCalibrator->UpdateCalHist(fCalMode);
+ }
}
+
// Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method
// DeleteHistos();
if (fdMemoryTime > 0.) CleanLHMemory();
@@ -621,7 +687,10 @@ Bool_t CbmTofEventClusterizer::RegisterInputs()
LOG(error) << GetName() << ": No Tof digi input!";
return kFALSE;
}
- if (fDigiMan->IsPresent(ECbmModuleId::kT0)) { LOG(info) << GetName() << ": found separate T0 digi input!"; }
+ if (fDigiMan->IsPresent(ECbmModuleId::kT0)) {
+ LOG(info) << GetName() << ": separate T0 digi input!";
+ //fT0DigiVec = fManager->InitObjectAs<std::vector<CbmTofDigi> const*>("TzdDigi");
+ }
else {
LOG(info) << "No separate T0 digi input found.";
} // if( ! fT0DigiVec )
@@ -632,10 +701,13 @@ Bool_t CbmTofEventClusterizer::RegisterInputs()
"TofTrbHeader Object";
}
- fEvtHeader = (FairEventHeader*) fManager->GetObject("EventHeader.");
- if (NULL == fEvtHeader) {
- LOG(info) << "CbmTofEventClusterizer::RegisterInputs => Could not get "
- "EvtHeader Object";
+ fEvtHeader = dynamic_cast<FairEventHeader*>(fManager->GetObject("EventHeader."));
+ if (NULL == fEvtHeader) { LOG(fatal) << "CbmTofEventClusterizer::RegisterInputs => Could not get EvtHeader Object"; }
+ fTsHeader = fManager->InitObjectAs<CbmTsEventHeader const*>("EventHeader."); //for data
+ if (NULL == fTsHeader) {
+ LOG(info) << "CbmTofEventClusterizer::RegisterInputs => Could not get TsHeader Object";
+ fTimeSlice = fManager->InitObjectAs<CbmTimeSlice const*>("TimeSlice."); // MC
+ if (NULL == fTimeSlice) { LOG(info) << "CbmTofEventClusterizer::RegisterInputs => Could not get TimeSlice Object"; }
}
if (NULL == fEventsColl) {
@@ -659,7 +731,7 @@ Bool_t CbmTofEventClusterizer::RegisterOutputs()
//fTofCalDigisColl = new TClonesArray("CbmTofDigi");
fTofCalDigiVec = new std::vector<CbmTofDigi>();
- fTofHitsColl = new TClonesArray("CbmTofHit");
+ fTofHitsColl = new TClonesArray("CbmTofHit", 100);
fTofDigiMatchColl = new TClonesArray("CbmMatch", 100);
@@ -685,7 +757,7 @@ Bool_t CbmTofEventClusterizer::RegisterOutputs()
fTofHitsCollOut = fTofHitsColl;
rootMgr->Register(tHitDigiMatchBranchName, "Tof", fTofDigiMatchColl, fbWriteHitsInOut);
- LOG(info) << Form("Register DigiMatch in TS mode at %p", fTofDigiMatchColl);
+ LOG(info) << Form("Register DigiMatch in TS mode at %p with %d", fTofDigiMatchColl, (int) fbWriteHitsInOut);
fTofDigiMatchCollOut = fTofDigiMatchColl;
if (NULL != fTofDigiPointMatches) {
@@ -703,7 +775,7 @@ Bool_t CbmTofEventClusterizer::RegisterOutputs()
rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVecOut, fbWriteDigisInOut);
rootMgr->Register(tHitBranchName, "Tof", fTofHitsCollOut, fbWriteHitsInOut);
rootMgr->Register(tHitDigiMatchBranchName, "Tof", fTofDigiMatchCollOut, fbWriteHitsInOut);
- LOG(info) << Form("Register DigiMatch in event mode at %p", fTofDigiMatchCollOut);
+ LOG(info) << Form("Register DigiMatch in event mode at %p with %d ", fTofDigiMatchCollOut, (int) fbWriteHitsInOut);
if (NULL != fTofDigiPointMatches) {
fTofDigiPointMatchesOut = new std::vector<CbmMatch>();
@@ -794,6 +866,9 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
fvCPTotOff.resize(iNbSmTypes);
fvCPWalk.resize(iNbSmTypes);
fvCPDelTof.resize(iNbSmTypes);
+ fvCPTOffY.resize(iNbSmTypes);
+ fvCPTOffYBinWidth.resize(iNbSmTypes);
+ fvCPTOffYRange.resize(iNbSmTypes);
for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
@@ -801,6 +876,9 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
fvCPTotOff[iSmType].resize(iNbSm * iNbRpc);
fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffY[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffYBinWidth[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffYRange[iSmType].resize(iNbSm * iNbRpc);
fvCPDelTof[iSmType].resize(iNbSm * iNbRpc);
for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
@@ -813,6 +891,8 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
for (Int_t iSel = 0; iSel < iNSel; iSel++)
fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] = 0.; // initialize
}
+ fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
+ fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
@@ -853,7 +933,7 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
// read parameter from histos
if (fCalParFileName.IsNull()) return kTRUE;
- fCalParFile = new TFile(fCalParFileName, "");
+ fCalParFile = new TFile(fCalParFileName, "READ");
if (NULL == fCalParFile) {
if (fCalMode % 10 == 9) { //modify reference file name
int iCalMode = fCalParFileName.Index("tofClust") - 3;
@@ -863,49 +943,45 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
if (NULL == fCalParFile)
LOG(fatal) << "CbmTofEventClusterizer::InitCalibParameter: "
<< "file " << fCalParFileName << " does not exist!";
+
return kTRUE;
}
+ LOG(fatal) << "Calibration parameter file not existing!";
}
- /*
- gDirectory->Print();
- fCalParFile->cd();
- fCalParFile->ls();
- */
+ fhSmCluSvel.resize(iNbSmTypes);
for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
Int_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
- TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
-
- // copy Histo to memory
- TDirectory* curdir = gDirectory;
- if (NULL != hSvel) {
- gDirectory->cd(oldDir->GetPath());
- // TProfile *hSvelmem = (TProfile *)hSvel->Clone(); (VF) not used
- gDirectory->cd(curdir->GetPath());
- }
- else {
- LOG(info) << "Svel histogram not found for module type " << iSmType;
- }
-
- for (Int_t iPar = 0; iPar < 4; iPar++) {
- TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Fpar%1d", iSmType, iPar));
- if (NULL != hFparcur) {
- gDirectory->cd(oldDir->GetPath());
- // TProfile *hFparmem = (TProfile *)hFparcur->Clone(); (VF) not used
- gDirectory->cd(curdir->GetPath());
+ if (iNbSm > 0 && NULL == fhSmCluSvel[iSmType]) {
+ TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
+ if (NULL == hSvel) LOG(warn) << "hSvel not found in " << gDirectory->GetName();
+ else {
+ fhSmCluSvel[iSmType] = (TProfile*) hSvel->Clone();
+ LOG(info) << fhSmCluSvel[iSmType]->GetName() << " cloned from " << gDirectory->GetName();
}
}
-
+ TDirectory* curdir = gDirectory;
for (Int_t iSm = 0; iSm < iNbSm; iSm++)
for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
-
// update default parameter
- if (NULL != hSvel) {
- Double_t Vscal = 1.; //hSvel->GetBinContent(iSm*iNbRpc+iRpc+1);
+ if (NULL != fhSmCluSvel[iSmType]) {
+ Double_t Vscal = fhSmCluSvel[iSmType]->GetBinContent(iSm * iNbRpc + iRpc + 1);
if (Vscal == 0.) Vscal = 1.;
+ if (Vscal < 0.8) {
+ LOG(warn) << "Fix parameter value of signal velocity TSR " << iSmType << iSm << iRpc << ": " << Vscal;
+ Vscal = 0.8;
+ }
+ if (Vscal > 1.2) {
+ LOG(warn) << "Fix parameter value of signal velocity TSR " << iSmType << iSm << iRpc << ": " << Vscal;
+ Vscal = 1.2;
+ }
+ Vscal *= fdModifySigvel; //1.03; // testing the effect of wrong signal velocity, FIXME
fDigiBdfPar->SetSigVel(iSmType, iSm, iRpc, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * Vscal);
- LOG(debug) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
- << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ LOG(info) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
+ << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ }
+ else {
+ LOG(warn) << "Svel modifier histo not found in " << curdir->GetName();
}
TH2D* htempPos_pfx =
(TH2D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
@@ -919,13 +995,11 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
//Int_t iNbinTot = htempTot_Mean->GetNbinsX();//not used any more
for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
-
for (Int_t iSide = 0; iSide < 2; iSide++) {
Double_t TotMean = htempTot_Mean->GetBinContent(iCh * 2 + 1 + iSide); //nh +1 empirical(?)
if (0.001 < TotMean) { fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] *= fdTTotMean / TotMean; }
fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
}
-
Double_t YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1); //nh +1 empirical(?)
Double_t TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
if (std::isnan(YMean) || std::isnan(TMean)) {
@@ -950,8 +1024,8 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
}
//if(iSmType==0 && iSm==4 && iRpc==2 && iCh==26)
- //if (iSmType==0 && iSm==1 && iRpc == 1)
- if (iSmType == 9 && iSm == 0 && iRpc == 0 && iCh == 10) // select specific channel
+ if (iSmType == 0 && iSm == 2 && iRpc == 0)
+ //if (iSmType == 9 && iSm == 0 && iRpc == 0 && iCh == 10) // select specific channel
LOG(info) << "InitCalibParameter:"
<< " TSRC " << iSmType << iSm << iRpc << iCh
<< Form(": YMean %6.3f, TYOff %6.3f, TMean %6.3f", YMean, dTYOff, TMean) << " -> "
@@ -965,6 +1039,20 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
+ if (NULL == htempWalk0 && NULL == htempWalk1) { // regenerate Walk histos
+ int iSide = 0;
+ htempWalk0 =
+ new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType,
+ iSm, iRpc, iCh, iSide),
+ nbClWalkBinX, fdTOTMin, fdTOTMax);
+ iSide = 1;
+ htempWalk1 =
+ new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType,
+ iSm, iRpc, iCh, iSide),
+ nbClWalkBinX, fdTOTMin, fdTOTMax);
+ }
if (NULL != htempWalk0 && NULL != htempWalk1) { // reinitialize Walk array
LOG(debug) << "Initialize Walk correction for "
<< Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
@@ -974,9 +1062,11 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
for (Int_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
- if (0) // iCh == 5 && iBin == 10) // debugging
- LOG(info) << Form("Read New SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d walk %f ", iSmType, iSm, iRpc, iCh,
- iBin, fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin]);
+ if (iSmType == 0 && iSm == 0 && iRpc == 2 && iCh == 15) // debugging
+ LOG(info) << Form("Read new SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d cen %f walk %f %f", iSmType, iSm,
+ iRpc, iCh, iBin, htempWalk0->GetBinCenter(iBin + 1),
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin],
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin]);
if (5 == iSmType || 8 == iSmType) { // Pad structure, enforce consistency
if (fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin]
!= fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin]) {
@@ -994,6 +1084,22 @@ Bool_t CbmTofEventClusterizer::InitCalibParameter()
LOG(info) << "No Walk histograms for TSRC " << iSmType << iSm << iRpc << iCh;
}
}
+ // look for TcorY corrections
+ LOG(info) << "Check for TCorY in " << gDirectory->GetName();
+ TH1* hTCorY =
+ (TH1*) gDirectory->FindObjectAny(Form("calXY_SmT%d_sm%03d_rpc%03d_TOff_z_all_TcorY", iSmType, iSm, iRpc));
+ if (NULL != hTCorY) {
+ fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = hTCorY->GetBinWidth(0);
+ fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = hTCorY->GetXaxis()->GetXmax();
+ LOG(info) << "Initialize TCorY: TSR " << iSmType << iSm << iRpc << ", Bwid "
+ << fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] << ", Range "
+ << fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc];
+ fvCPTOffY[iSmType][iSm * iNbRpc + iRpc].resize(hTCorY->GetNbinsX());
+ for (int iB = 0; iB < hTCorY->GetNbinsX(); iB++) {
+ //LOG(info) << "iB " << iB <<": " << hTCorY->GetBinContent(iB+1);
+ fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][iB] = hTCorY->GetBinContent(iB + 1);
+ }
+ }
}
else {
LOG(warning) << " Calibration histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
@@ -1038,7 +1144,7 @@ Bool_t CbmTofEventClusterizer::LoadGeometry()
{
LOG(info) << "CbmTofEventClusterizer::LoadGeometry starting for " << fDigiBdfPar->GetNbDet()
<< " described detectors, " << fDigiPar->GetNrOfModules() << " geometrically known cells ";
-
+ fCellIdGeoMap.clear();
Int_t iNrOfCells = fDigiPar->GetNrOfModules();
LOG(info) << "Digi Parameter container contains " << iNrOfCells << " cells.";
for (Int_t icell = 0; icell < iNrOfCells; ++icell) {
@@ -1064,6 +1170,8 @@ Bool_t CbmTofEventClusterizer::LoadGeometry()
gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
LOG(debug2) << Form(" Node at (%6.1f,%6.1f,%6.1f) : 0x%p", fChannelInfo->GetX(), fChannelInfo->GetY(),
fChannelInfo->GetZ(), fNode);
+ fCellIdGeoMap[cellId] = gGeoManager->MakePhysicalNode();
+
if (icell == 0) {
TGeoHMatrix* cMatrix = gGeoManager->GetCurrentMatrix();
fNode->Print();
@@ -1083,8 +1191,8 @@ Bool_t CbmTofEventClusterizer::LoadGeometry()
Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
Int_t iSmId = CbmTofAddress::GetSmId(iUniqueId);
Int_t iRpcId = CbmTofAddress::GetRpcId(iUniqueId);
- LOG(info) << " DetIndx " << iDetIndx << "(" << iNbDet << "), SmType " << iSmType << ", SmId " << iSmId << ", RpcId "
- << iRpcId << " => UniqueId " << Form("0x%08x ", iUniqueId)
+ LOG(info) << " DetIndx " << iDetIndx << "(" << iNbDet << "), TSR " << iSmType << iSmId << iRpcId << " => UniqueId "
+ << Form("0x%08x ", iUniqueId)
<< Form(" Svel %6.6f, DeadStrips 0x%08x ", fDigiBdfPar->GetSigVel(iSmType, iSmId, iRpcId),
fvDeadStrips[iDetIndx]);
Int_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpcId);
@@ -1195,6 +1303,10 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
fhClustBuildTime =
new TH1I("TofEventClus_ClustBuildTime", "Time needed to build clusters in each event; Time [s]", 4000, 0.0, 4.0);
+
+ fhClustHitsDigi =
+ new TH2D(Form("hClustHitsDigi"), Form("Hits vs. CalDigis; Mul_{Digi}; Mul_{Hits}"), 100, 0, 500, 100, 0, 200);
+
// RPC related distributions
Int_t iNbDet = fDigiBdfPar->GetNbDet();
fviDetId.resize(iNbDet);
@@ -1206,7 +1318,14 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
fviDetId[iDetIndx] = iUniqueId;
}
- if (fDutId < 0) return kTRUE;
+ if (fTotMax != 0.) fdTOTMax = fTotMax;
+ if (fTotMin != 0.) fdTOTMin = fTotMin;
+ LOG(info) << "ToT init to Min " << fdTOTMin << " Max " << fdTOTMax;
+
+ if (fDutId < 0) {
+ LOG(info) << GetName() << ": Skip booking of calibration histograms ";
+ return kTRUE;
+ }
// Sm related distributions
fhSmCluPosition.resize(fDigiBdfPar->GetNbSmTypes());
@@ -1258,16 +1377,21 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 99, -TSumMax * 2., TSumMax * 2.);
- TProfile* hSvelcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iS));
- if (NULL == hSvelcur) {
- fhSmCluSvel[iS] =
- new TProfile(Form("cl_SmT%01d_Svel", iS), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iS),
- fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
- fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0.8, 1.2);
- fhSmCluSvel[iS]->Sumw2();
+ if (NULL == fhSmCluSvel[iS]) {
+ TProfile* hSvelcur = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iS));
+ if (NULL == hSvelcur) {
+ LOG(info) << "Histo " << Form("cl_SmT%01d_Svel", iS) << " not found, recreate ...";
+ fhSmCluSvel[iS] =
+ new TProfile(Form("cl_SmT%01d_Svel", iS), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iS),
+ fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0,
+ fDigiBdfPar->GetNbSm(iS) * fDigiBdfPar->GetNbRpc(iS), 0.8, 1.2);
+ fhSmCluSvel[iS]->Sumw2();
+ }
+ else {
+ fhSmCluSvel[iS] = (TProfile*) hSvelcur->Clone();
+ LOG(info) << fhSmCluSvel[iS]->GetName() << " cloned from " << gDirectory->GetName();
+ }
}
- else
- fhSmCluSvel[iS] = (TProfile*) hSvelcur->Clone();
fhSmCluFpar[iS].resize(4);
for (Int_t iPar = 0; iPar < 4; iPar++) {
@@ -1316,6 +1440,10 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
fhRpcDigiDTFD.resize(iNbDet);
fhRpcDigiDTMul.resize(iNbDet);
fhRpcDigiRate.resize(iNbDet);
+ fhRpcDigiTotLeft.resize(iNbDet);
+ fhRpcDigiTotRight.resize(iNbDet);
+ fhRpcDigiTotDiff.resize(iNbDet);
+ fhRpcDigiTotMap.resize(iNbDet);
fhRpcCluMul.resize(iNbDet);
fhRpcCluRate.resize(iNbDet);
fhRpcCluRate10s.resize(iNbDet);
@@ -1338,9 +1466,6 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
fhRpcDTLastHits_Tot.resize(iNbDet);
fhRpcDTLastHits_CluSize.resize(iNbDet);
- if (fTotMax != 0.) fdTOTMax = fTotMax;
- if (fTotMin != 0.) fdTOTMin = fTotMin;
-
for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
@@ -1459,6 +1584,25 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
Form("cl_SmT%01d_sm%03d_rpc%03d_GloPos", iSmType, iSmId, iRpcId),
Form("Clu global position of Rpc #%03d in Sm %03d of type %d; xpos [cm]; ypos [cm]", iRpcId, iSmId, iSmType), 200,
-100, 100, 200, -100, 100);
+ fhRpcDigiTotLeft[iDetIndx] =
+ new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_Left", iSmType, iSmId, iRpcId),
+ Form("Digi Tot of Rpc #%03d in Sm %03d of type %d; Tot; ypos [cm]", iRpcId, iSmId, iSmType), 100,
+ fdTOTMin, fdTOTMax, 99, -YDMAX, YDMAX);
+
+ fhRpcDigiTotRight[iDetIndx] =
+ new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_Right", iSmType, iSmId, iRpcId),
+ Form("Digi Tot of Rpc #%03d in Sm %03d of type %d; Tot; ypos [cm]", iRpcId, iSmId, iSmType), 100,
+ fdTOTMin, fdTOTMax, 99, -YDMAX, YDMAX);
+
+ fhRpcDigiTotDiff[iDetIndx] =
+ new TH2D(Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_Diff", iSmType, iSmId, iRpcId),
+ Form("Digi Tot of Rpc #%03d in Sm %03d of type %d; Tot; ypos [cm]", iRpcId, iSmId, iSmType), 200,
+ fdTOTMax, fdTOTMax, 99, -YDMAX, YDMAX);
+
+ fhRpcDigiTotMap[iDetIndx] = new TH2D(
+ Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_Map", iSmType, iSmId, iRpcId),
+ Form("Digi Tot left vs Tot right of Rpc #%03d in Sm %03d of type %d; Tot; ypos [cm]", iRpcId, iSmId, iSmType),
+ 100, 0, fdTOTMax, 100, 0, fdTOTMax);
fhRpcCluPositionEvol[iDetIndx] = new TProfile(Form("cl_SmT%01d_sm%03d_rpc%03d_PosEvol", iSmType, iSmId, iRpcId),
Form("Clu position of Rpc #%03d in Sm %03d of type %d; Analysis Time "
@@ -1830,6 +1974,7 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
"from a single TofPoint; tMcPoint -tTofHit [ps]",
2000, 0.0, 50.0, 2000, 0.0, 50.0);
} // else of if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
+ /*
fhClusterSize = new TH1I("Clus_ClusterSize", "Cluster Size distribution; Cluster Size [Strips]", 100, 0.5, 100.5);
fhClusterSizeType =
new TH2I("Clus_ClusterSizeType",
@@ -1856,9 +2001,9 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
fhMultiTrkProbPos = new TH2D("Clus_MultiTrkProbPos",
"Probability of having a cluster with multiple tracks as "
"function of position; X [cm]; Y [cm]; Prob. [%]",
- 1500, -750, 750, 1000, -500, 500);
+ 1500, -750, 750, 1000, -500, 500);
} // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
-
+ */
gDirectory->cd(oldir->GetPath()); // <= To prevent histos from being sucked in by the param file of the TRootManager!
return kTRUE;
@@ -1867,6 +2012,7 @@ Bool_t CbmTofEventClusterizer::CreateHistos()
Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
{
fhClustBuildTime->Fill(fStop.GetSec() - fStart.GetSec() + (fStop.GetNanoSec() - fStart.GetNanoSec()) / 1e9);
+ fhClustHitsDigi->Fill(fTofCalDigiVec->size(), fTofHitsColl->GetEntriesFast());
if (fDutId < 0) return kTRUE;
@@ -1980,7 +2126,8 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
//LOG(debug)<<"TimeAna "<<StartAnalysisTime<<", "<< pHit->GetTime()<<", "<<dTimeAna;
fhRpcCluRate[iDetIndx]->Fill(dTimeAna, 1. / fhRpcCluRate[iDetIndx]->GetBinWidth(1));
- // deal with spill structure
+ // deal with spill structures
+
Double_t dTimeAna10s = pHit->GetTime() + dTsStartTime - fdStartAna10s;
if (iHitInd == 0) {
fhEvCluMul->Fill(dTimeAna, dCluMul);
@@ -2337,7 +2484,6 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
uTriggerPattern |= (0x1 << (iSel * 3 + CbmTofAddress::GetRpcId(pTrig[iSel]->GetAddress() & DetMask)));
}
}
-
LOG(debug1) << "Inspect trigger pattern";
for (Int_t iSel = 0; iSel < iNSel; iSel++) {
if (BSel[iSel]) {
@@ -2389,7 +2535,7 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
continue;
}
/*TGeoNode *fNode=*/ // prepare global->local trafo
- gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
LOG(debug1) << "Hit info: "
<< Form(" 0x%08x %d %f %f %f %f %f %d", iChId, iCh, pHit->GetX(), pHit->GetY(), pHit->GetTime(),
@@ -2400,12 +2546,14 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
hitpos[1] = pHit->GetY();
hitpos[2] = pHit->GetZ();
Double_t hitpos_local[3];
- TGeoNode* cNode = gGeoManager->GetCurrentNode();
- gGeoManager->MasterToLocal(hitpos, hitpos_local);
+ //TGeoNode* cNode = gGeoManager->GetCurrentNode();
+ fCellIdGeoMap[iChId]->GetMatrix()->MasterToLocal(hitpos, hitpos_local);
+ /*
LOG(debug1) << Form(" MasterToLocal for %d, %d%d%d, node %p: "
"(%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
iDetIndx, iSmType, iSm, iRpc, cNode, hitpos[0], hitpos[1], hitpos[2], hitpos_local[0],
hitpos_local[1], hitpos_local[2]);
+ */
Bool_t bFillPos = kTRUE;
//if( fCalMode/10 > 4 && pHit->GetClusterSize() < 3 ) bFillPos=kFALSE;
if (bFillPos) {
@@ -2511,6 +2659,11 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
CbmTofDigi* pDigS1 = pDig1;
if (iS1 == 0) pDigS1 = pDig0;
+ fhRpcDigiTotLeft[iDetIndx]->Fill(pDigS0->GetTot(), hitpos_local[1]);
+ fhRpcDigiTotRight[iDetIndx]->Fill(pDigS1->GetTot(), hitpos_local[1]);
+ fhRpcDigiTotDiff[iDetIndx]->Fill(pDigS0->GetTot() - pDigS1->GetTot(), hitpos_local[1]);
+ fhRpcDigiTotMap[iDetIndx]->Fill(pDigS0->GetTot(), pDigS1->GetTot());
+
if (iCh0 != iCh1 || iS0 == iS1) {
LOG(error) << Form(" MT2 for Tofhit %d in iDetIndx %d, Ch %d from "
"in event %d, ",
@@ -2834,46 +2987,47 @@ Bool_t CbmTofEventClusterizer::FillHistos(CbmEvent* tEvent)
}
} // iHitInd hit loop end
- for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
- for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
- LOG(debug1) << "CbmTofEventClusterizer::FillHistos: "
- << Form(" %3d %3d %3lu ", iSmType, iRpc, fviClusterSize[iSmType][iRpc].size());
+ if (false)
+ for (Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) {
+ for (Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
+ LOG(debug1) << "CbmTofEventClusterizer::FillHistos: "
+ << Form(" %3d %3d %3lu ", iSmType, iRpc, fviClusterSize[iSmType][iRpc].size());
- for (UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++) {
- LOG(debug2) << "CbmTofEventClusterizer::FillHistos: " << Form(" %3d %3d %3d ", iSmType, iRpc, uCluster);
+ for (UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++) {
+ LOG(debug2) << "CbmTofEventClusterizer::FillHistos: " << Form(" %3d %3d %3d ", iSmType, iRpc, uCluster);
- fhClusterSize->Fill(fviClusterSize[iSmType][iRpc][uCluster]);
- fhClusterSizeType->Fill(fviClusterSize[iSmType][iRpc][uCluster],
- 40 * iSmType + iRpc); //FIXME - hardwired constant
- if (kFALSE) // kTRUE == fDigiBdfPar->ClustUseTrackId() )
- {
- fhTrackMul->Fill(fviTrkMul[iSmType][iRpc][uCluster]);
- fhClusterSizeMulti->Fill(fviClusterSize[iSmType][iRpc][uCluster], fviTrkMul[iSmType][iRpc][uCluster]);
- if (1 == fviTrkMul[iSmType][iRpc][uCluster])
- fhTrk1MulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
- if (1 < fviTrkMul[iSmType][iRpc][uCluster])
- fhHiTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
- fhAllTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
- } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
- if (kFALSE) // 1 == fviTrkMul[iSmType][iRpc][uCluster] )
- {
- fhClustSizeDifX->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
- fhClustSizeDifY->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
- if (1 == fviClusterSize[iSmType][iRpc][uCluster]) {
- fhChDifDifX->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
- fhChDifDifY->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
+ fhClusterSize->Fill(fviClusterSize[iSmType][iRpc][uCluster]);
+ fhClusterSizeType->Fill(fviClusterSize[iSmType][iRpc][uCluster],
+ 40 * iSmType + iRpc); //FIXME - hardwired constant
+ if (kFALSE) // kTRUE == fDigiBdfPar->ClustUseTrackId() )
+ {
+ fhTrackMul->Fill(fviTrkMul[iSmType][iRpc][uCluster]);
+ fhClusterSizeMulti->Fill(fviClusterSize[iSmType][iRpc][uCluster], fviTrkMul[iSmType][iRpc][uCluster]);
+ if (1 == fviTrkMul[iSmType][iRpc][uCluster])
+ fhTrk1MulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
+ if (1 < fviTrkMul[iSmType][iRpc][uCluster])
+ fhHiTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
+ fhAllTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster]);
+ } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() )
+ if (kFALSE) // 1 == fviTrkMul[iSmType][iRpc][uCluster] )
+ {
+ fhClustSizeDifX->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
+ fhClustSizeDifY->Fill(fviClusterSize[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
+ if (1 == fviClusterSize[iSmType][iRpc][uCluster]) {
+ fhChDifDifX->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]);
+ fhChDifDifY->Fill(fvdDifCh[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]);
+ }
}
- }
- } // for( UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++ )
- fviClusterSize[iSmType][iRpc].clear();
- fviTrkMul[iSmType][iRpc].clear();
- fvdX[iSmType][iRpc].clear();
- fvdY[iSmType][iRpc].clear();
- fvdDifX[iSmType][iRpc].clear();
- fvdDifY[iSmType][iRpc].clear();
- fvdDifCh[iSmType][iRpc].clear();
- } // for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ )
- }
+ } // for( UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++ )
+ fviClusterSize[iSmType][iRpc].clear();
+ fviTrkMul[iSmType][iRpc].clear();
+ fvdX[iSmType][iRpc].clear();
+ fvdY[iSmType][iRpc].clear();
+ fvdDifX[iSmType][iRpc].clear();
+ fvdDifY[iSmType][iRpc].clear();
+ fvdDifCh[iSmType][iRpc].clear();
+ } // for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ )
+ }
fhNbSameSide->Fill(fiNbSameSide);
} // if(0<iNbTofHits) end
@@ -3777,16 +3931,16 @@ Bool_t CbmTofEventClusterizer::WriteHistos()
fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
}
/*
- Double_t TotMean=((TProfile *)htempTot_pfx)->GetBinContent(iCh+1);
- if(0.001 < TotMean){
- fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0] *= fdTTotMean / TotMean;
- fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1] *= fdTTotMean / TotMean;
- }
- */
+ Double_t TotMean=((TProfile *)htempTot_pfx)->GetBinContent(iCh+1);
+ if(0.001 < TotMean){
+ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0] *= fdTTotMean / TotMean;
+ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1] *= fdTTotMean / TotMean;
+ }
+ */
if (fCalMode < 90) // keep digi TOT calibration in last step
for (Int_t iSide = 0; iSide < 2; iSide++) {
Int_t ib = iCh * 2 + 1 + iSide;
- TH1* hbin = htempTot->ProjectionY(Form("bin%d", ib), ib, ib);
+ TH1* hbin = htempTot->ProjectionY(Form("%s_bin%d", htempTot->GetName(), ib), ib, ib);
if (100 > hbin->GetEntries()) continue; //request min number of entries
/* Double_t Ymax=hbin->GetMaximum();*/
Int_t iBmax = hbin->GetMaximumBin();
@@ -3800,6 +3954,38 @@ Bool_t CbmTofEventClusterizer::WriteHistos()
}
// Double_t TotMean=htempTot_Mean->GetBinContent(ib);
Double_t TotMean = hbin->GetMean();
+ Double_t dCtsTot = hbin->GetEntries();
+ if (dCtsTot > WalkNHmin) {
+ if (iBmax == 0) {
+ if (hbin->GetBinContent(hbin->GetNbinsX() > 0))
+ TotMean = xaxis->GetBinCenter(hbin->GetNbinsX() - 1);
+ LOG(info) << "TotInvalid " << hbin->GetName() << ": " << TotMean;
+ }
+ double dFMean = hbin->GetBinCenter(iBmax);
+ double dFLim = dFMean * 0.5; // CAUTION, hardcoded parameter
+ LOG(info) << Form("FitTot TSRC %d%d%d%2d: Mean %6.2f, Width %6.2f ", iSmType, iSm, iRpc, ib,
+ dFMean, dFLim);
+ if (dFMean > 2.) {
+ TFitResultPtr fRes = hbin->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ TotMean = fRes->Parameter(1); //overwrite mean
+ }
+ else {
+ TotMean = dFMean;
+ LOG(warn) << "TotFitFail " << hbin->GetName() << ": " << TotMean << ", wid " << dFLim << ", "
+ << gMinuit->fCstatu;
+ }
+ }
+ else {
+ TotMean = dFMean;
+ if (TotMean < 0.2)
+ //if(hbin->GetBinContent(hbin->GetNbinsX()) > 0 )
+ TotMean = xaxis->GetBinCenter(hbin->GetNbinsX() - 1);
+ LOG(warn) << "TotooR " << hbin->GetName() << ": " << TotMean << ", gain "
+ << fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] << ", ovfl "
+ << hbin->GetBinContent(hbin->GetNbinsX());
+ }
+ }
if (15 == iSmType) {
LOG(warning) << "Gain for "
<< Form("SmT%01d_sm%03d_rpc%03d_Side%d: TotMean %f, prof %f, "
@@ -3824,13 +4010,17 @@ Bool_t CbmTofEventClusterizer::WriteHistos()
Double_t dTOff =
0.5
* (fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] + fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
+ /*
Double_t dGain = 0.5
* (fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0]
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1]);
+ */
+ Double_t dGain0 = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0];
+ Double_t dGain1 = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1];
fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] = dTOff;
fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] = dTOff;
- fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0] = dGain;
- fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] = dGain;
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][0] = dGain0;
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] = dGain1;
} // diamond warning end
} // for( Int_t iCh = 0; iCh < iNbCh; iCh++ )
} // iSmType selection condition
@@ -4080,7 +4270,7 @@ Bool_t CbmTofEventClusterizer::WriteHistos()
{
Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
- Int_t iYCalMode = 1;
+ Int_t iYCalMode = 0; // FIXME hardcoded option
if ((fCalSmAddr < 0) || (fCalSmAddr != iSmAddr)) { // select detectors for updating offsets
LOG(debug) << "WriteHistos (calMode==5): update Offsets and Gains, "
"keep Walk and DelTof for "
@@ -4543,7 +4733,7 @@ Bool_t CbmTofEventClusterizer::WriteHistos()
// fhTSmCluTOff[iS][iSel]->Write();
// fhTSmCluTRun[iS][iSel]->Write();
}
- } // iFindT0 loop end
+ }
/// Restore old global file and folder pointer to avoid messing with FairRoot
gFile = oldFile;
@@ -4600,6 +4790,7 @@ Bool_t CbmTofEventClusterizer::BuildClusters()
return kFALSE;
}
fiNevtBuild++;
+
LOG(debug) << "Build clusters from "
// <<fTofDigisColl->GetEntriesFast()<<" digis in event "<<fiNevtBuild;
<< fTofDigiVec.size() << " digis in event " << fiNevtBuild;
@@ -4841,8 +5032,7 @@ Bool_t CbmTofEventClusterizer::BuildClusters()
// Then loop over the digis array and store the Digis in separate vectors for
// each RPC modules
-
- // iNbTofDigi = fTofCalDigisColl->GetEntriesFast();
+ // iNbTofDigi = fTofCalDigisColl->GetEntriesFast();
iNbTofDigi = fTofCalDigiVec->size();
for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
// pDigi = (CbmTofDigi*) fTofCalDigisColl->At( iDigInd );
@@ -4915,8 +5105,9 @@ Bool_t CbmTofEventClusterizer::BuildClusters()
}
}
}
- BuildHits();
+ BuildHits();
+ CalibHits();
} // if( kTRUE ) obsolete )
return kTRUE;
@@ -5234,8 +5425,21 @@ Bool_t CbmTofEventClusterizer::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc,
CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
iChId = fTofId->SetDetectorInfo(xDetInfo);
fChannelInfo = fDigiPar->GetCell(iChId);
- gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
-
+ if (iChId == 0x13800036) {
+ const double local[3] = {3 * 0};
+ double global[3];
+ fCellIdGeoMap[iChId]->GetMatrix()->LocalToMaster(local, global);
+ LOG(info) << "GeoTest1: " << global[0] << ", " << global[1] << ", " << global[2];
+ }
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ()); //needed, why?
+ if (iChId == 0x13800036) {
+ TGeoHMatrix* pMatrix = gGeoManager->GetCurrentMatrix();
+ const double local[3] = {3 * 0};
+ double global[3];
+ pMatrix->LocalToMaster(local, global);
+ LOG(info) << "GeoTest2: " << global[0] << ", " << global[1] << ", " << global[2];
+ if (fiNbHits == 10) LOG(fatal) << " for inspection";
+ }
Double_t dTimeDif = xDigiA->GetTime() - xDigiB->GetTime();
Double_t dPosY = 0.;
if (1 == xDigiA->GetSide()) dPosY = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * dTimeDif * 0.5;
@@ -5302,18 +5506,21 @@ Bool_t CbmTofEventClusterizer::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc,
}
*/
Double_t hitpos[3] = {3 * 0.};
- if (5 != iSmType) { // Diamond beam counter always at (0,0,0)
- /*TGeoNode* cNode = */ gGeoManager->GetCurrentNode();
- /*TGeoHMatrix* cMatrix = */ gGeoManager->GetCurrentMatrix();
- gGeoManager->LocalToMaster(hitpos_local, hitpos);
- }
- TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
- TVector3 hitPosErr(0.5, 0.5, 0.5); // FIXME including positioning uncertainty
Int_t iChm = floor(dLastPosX / fChannelInfo->GetSizex()) + iNbCh / 2;
if (iChm < 0) iChm = 0;
if (iChm > iNbCh - 1) iChm = iNbCh - 1;
iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iChm, 0, iSmType);
+ if (5 != iSmType) { // Diamond beam counter always at (0,0,0)
+ fCellIdGeoMap[iDetId]->GetMatrix()->LocalToMaster(hitpos_local, hitpos);
+ /*TGeoNode* cNode = */ // gGeoManager->GetCurrentNode();
+ /*TGeoHMatrix* cMatrix = */ // gGeoManager->GetCurrentMatrix();
+ //gGeoManager->LocalToMaster(hitpos_local, hitpos);
+ }
+ TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
+ TVector3 hitPosErr(0.5, 0.5, 0.5); // FIXME including positioning uncertainty
+
+
Int_t iNbChanInHit = vDigiIndRef.size() / 2;
TString cstr = "Save A-Hit ";
@@ -5483,7 +5690,7 @@ Bool_t CbmTofEventClusterizer::BuildHits()
<< Form(" %3d %3d %3d %3d %3d ", iSmType, iSm, iRpc, fDigiBdfPar->GetChanOrient(iSmType, iRpc),
iNbCh);
- if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
+ if (0) { //1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) { // option disabled, used in Calibrator
// Horizontal strips => X comes from left right time difference
} // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
else { // Vertical strips => Y comes from bottom top time difference
@@ -5569,27 +5776,34 @@ Bool_t CbmTofEventClusterizer::BuildHits()
/* Int_t iLastChId = iChId; // Save Last hit channel*/
// 2 Digis = both sides present
- CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
- iChId = fTofId->SetDetectorInfo(xDetInfo);
- Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
+ //CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
+ //iChId = fTofId->SetDetectorInfo(xDetInfo);
+ iChId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
+ /*
LOG(debug1) << Form(" TSRC %d%d%d%d size %3lu ", iSmType, iSm, iRpc, iCh,
fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh].size())
<< Form(" ChId: 0x%08x 0x%08x ", iChId, iUCellId);
+ */
fChannelInfo = fDigiPar->GetCell(iChId);
if (NULL == fChannelInfo) {
- LOG(error) << "CbmTofEventClusterizer::BuildClusters: no "
- "geometry info! "
- << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ", iSmType, iSm, iRpc, iCh, iChId, iUCellId);
+ LOG(error) << "CbmTofEventClusterizer::BuildClusters: no geometry info! "
+ << Form(" %3d %3d %3d %3d 0x%08x", iSmType, iSm, iRpc, iCh, iChId);
break;
}
- TGeoNode* fNode = // prepare local->global trafo
- gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ //TGeoNode* fNode = // prepare local->global trafo
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ()); //needed for position spectra, why?
+ if (iChId == 0x13800036) {
+ const double local[3] = {3 * 0};
+ double global[3];
+ fCellIdGeoMap[iChId]->GetMatrix()->LocalToMaster(local, global);
+ LOG(info) << "GeoTest3: " << global[0] << ", " << global[1] << ", " << global[2];
+ } /*
LOG(debug2) << Form(" Node at (%6.1f,%6.1f,%6.1f) : %p", fChannelInfo->GetX(), fChannelInfo->GetY(),
fChannelInfo->GetZ(), fNode);
// fNode->Print();
-
+ */
CbmTofDigi* xDigiA = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][0];
CbmTofDigi* xDigiB = fStorDigi[iSmType][iSm * iNbRpc + iRpc][iCh][1];
@@ -5746,12 +5960,12 @@ Bool_t CbmTofEventClusterizer::BuildHits()
*/
Double_t hitpos[3] = {3 * 0.};
if (5 != iSmType) {
- /*TGeoNode* cNode =*/gGeoManager->GetCurrentNode();
- /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
+ fCellIdGeoMap[iChId]->GetMatrix()->LocalToMaster(hitpos_local, hitpos);
+ /*TGeoNode* cNode =*/ //gGeoManager->GetCurrentNode();
+ /*TGeoHMatrix* cMatrix =*/ //gGeoManager->GetCurrentMatrix();
//cNode->Print();
//cMatrix->Print();
-
- gGeoManager->LocalToMaster(hitpos_local, hitpos);
+ //gGeoManager->LocalToMaster(hitpos_local, hitpos);
}
LOG(debug1) << Form(" LocalToMaster: (%6.1f,%6.1f,%6.1f) "
"->(%6.1f,%6.1f,%6.1f)",
@@ -5879,12 +6093,13 @@ Bool_t CbmTofEventClusterizer::BuildHits()
}
fiNbHits++;
- // For Histogramming
+ // For Histogramming, vectors are cleared in FillHistos!
+ /*
fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
//fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
fvdX[iSmType][iRpc].push_back(dWeightedPosX);
fvdY[iSmType][iRpc].push_back(dWeightedPosY);
- /* no TofPoint available for data!
+ // no TofPoint available for data!
fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
@@ -6000,7 +6215,7 @@ Bool_t CbmTofEventClusterizer::BuildHits()
LOG(debug1) << "Process cluster " << iNbChanInHit;
// Check orientation to properly assign errors
- if (1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) {
+ if (0) { // 1 == fDigiBdfPar->GetChanOrient(iSmType, iRpc)) { //option disabled, LRF fixed
LOG(debug1) << "H-Hit ";
} // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) )
else {
@@ -6024,11 +6239,12 @@ Bool_t CbmTofEventClusterizer::BuildHits()
*/
Double_t hitpos[3] = {3 * 0.};
if (5 != iSmType) {
- /*TGeoNode* cNode=*/gGeoManager->GetCurrentNode();
- /*TGeoHMatrix* cMatrix =*/gGeoManager->GetCurrentMatrix();
+ fCellIdGeoMap[iChId]->GetMatrix()->LocalToMaster(hitpos_local, hitpos);
+ /*TGeoNode* cNode=*/ //gGeoManager->GetCurrentNode();
+ /*TGeoHMatrix* cMatrix =*/ //gGeoManager->GetCurrentMatrix();
//cNode->Print();
//cMatrix->Print();
- gGeoManager->LocalToMaster(hitpos_local, hitpos);
+ //gGeoManager->LocalToMaster(hitpos_local, hitpos);
}
LOG(debug1) << Form(" LocalToMaster for V-node: "
"(%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)",
@@ -6111,11 +6327,11 @@ Bool_t CbmTofEventClusterizer::BuildHits()
fiNbHits++;
// For Histogramming
+ /*
fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit);
//fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() );
fvdX[iSmType][iRpc].push_back(dWeightedPosX);
fvdY[iSmType][iRpc].push_back(dWeightedPosY);
- /*
fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX);
fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY);
fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan );
@@ -6220,6 +6436,7 @@ Bool_t CbmTofEventClusterizer::CalibRawDigis()
// walk correction
Double_t dTotBinSize = (fdTOTMax - fdTOTMin) / nbClWalkBinX;
Int_t iWx = (Int_t)((pCalDigi->GetTot() - fdTOTMin) / dTotBinSize);
+ // LOG(info)<<"Wx: "<<pCalDigi->GetTot()<<" "<<fdTOTMin<<" "<<fdTOTMax<<" "<<dTotBinSize<<" "<<iWx;
if (0 > iWx) iWx = 0;
if (iWx >= nbClWalkBinX) iWx = nbClWalkBinX - 1;
Double_t dDTot = (pCalDigi->GetTot() - fdTOTMin) / dTotBinSize - (Double_t) iWx - 0.5;
@@ -6254,9 +6471,10 @@ Bool_t CbmTofEventClusterizer::CalibRawDigis()
pCalDigi->SetTime(pCalDigi->GetTime() - dWT); // calibrate Digi Time
LOG(debug2) << " CluCal-Walk: " << pCalDigi->ToString();
- if (0) { //pDigi->GetType()==8) { // && pDigi->GetSm()==0){
+ if (
+ 0) { //pCalDigi->GetType()==0 && pCalDigi->GetSm()==0 && pCalDigi->GetRpc()==2 && pCalDigi->GetChannel()==15 ){
LOG(info)
- << "BuildClusters: CalDigi " << Form("%02d TSRCS ", iDigIndCal) << pCalDigi->GetType() << pCalDigi->GetSm()
+ << "CalDigi " << Form("%02d TSRCS ", iDigIndCal) << pCalDigi->GetType() << pCalDigi->GetSm()
<< pCalDigi->GetRpc() << Form("%02d", Int_t(pCalDigi->GetChannel())) << pCalDigi->GetSide()
<< Form(", T %15.3f", pCalDigi->GetTime()) << ", Tot " << pCalDigi->GetTot() << ", TotGain "
<< fvCPTotGain[pCalDigi->GetType()][pCalDigi->GetSm() * fDigiBdfPar->GetNbRpc(pCalDigi->GetType())
@@ -6368,7 +6586,7 @@ Bool_t CbmTofEventClusterizer::CalibRawDigis()
return kTRUE;
}
-void CbmTofEventClusterizer::SetDeadStrips(Int_t iDet, Int_t ival)
+void CbmTofEventClusterizer::SetDeadStrips(Int_t iDet, UInt_t ival)
{
if (fvDeadStrips.size() < static_cast<size_t>(iDet + 1)) fvDeadStrips.resize(iDet + 1);
fvDeadStrips[iDet] = ival;
@@ -6395,10 +6613,13 @@ double* CbmTofEventClusterizer::find_yedges(const char* hname, double dThr, doub
&& h1->GetBinContent(iBin + 2) > dLev) {
iBl = iBin;
xLow = h1->GetBinCenter(iBin);
+ break;
}
- if (iBl > -1) {
- if (iBh < 0 && h1->GetBinContent(iBin) < dLev && h1->GetBinContent(iBin + 1) < dLev
- && h1->GetBinContent(iBin + 2) < dLev) {
+ }
+ if (iBl > -1) {
+ for (int iBin = h1->GetNbinsX() - 1; iBin > 3; iBin--) {
+ if (iBh < 0 && h1->GetBinContent(iBin) > dLev && h1->GetBinContent(iBin - 1) > dLev
+ && h1->GetBinContent(iBin - 2) > dLev) {
iBh = iBin;
xHigh = h1->GetBinCenter(iBin);
break;
@@ -6413,6 +6634,7 @@ double* CbmTofEventClusterizer::find_yedges(const char* hname, double dThr, doub
dValid++;
}
} // for loop end
+ if (dValid == 0) { LOG(warn) << " No valid yEdge with thr " << dThr << " for " << h1->GetName(); }
dRes[0] = dLength;
dRes[1] = dMean;
}
@@ -6433,5 +6655,104 @@ CbmTofHit* CbmTofEventClusterizer::GetHitPointer(int iHitInd)
CbmMatch* CbmTofEventClusterizer::GetMatchPointer(int iHitInd)
{
+ if (fiHitStart + iHitInd > fTofDigiMatchCollOut->GetEntriesFast()) {
+ LOG(fatal) << "Invalid Hit index " << iHitInd << ", " << fiHitStart << ", "
+ << fTofDigiMatchCollOut->GetEntriesFast();
+ }
return (CbmMatch*) fTofDigiMatchCollOut->At(fiHitStart + iHitInd);
}
+
+CbmMatch* CbmTofEventClusterizer::GetMatchIndexPointer(int idx)
+{
+ if (idx > fTofDigiMatchCollOut->GetEntriesFast()) {
+ LOG(fatal) << "Invalid Hit index " << idx << ", " << fTofDigiMatchCollOut->GetEntriesFast();
+ }
+ return (CbmMatch*) fTofDigiMatchCollOut->At(idx);
+}
+
+double CbmTofEventClusterizer::GetLocalX(CbmTofHit* pHit)
+{
+ Int_t iChId = pHit->GetAddress();
+ Double_t hitpos[3];
+ Double_t hitpos_local[3] = {3 * 0.};
+ hitpos[0] = pHit->GetX();
+ hitpos[1] = pHit->GetY();
+ hitpos[2] = pHit->GetZ();
+ fChannelInfo = fDigiPar->GetCell(iChId);
+ Int_t iCh = CbmTofAddress::GetChannelId(iChId);
+ if (NULL == fChannelInfo) {
+ LOG(fatal) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
+ }
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ /*
+ gGeoManager->MasterToLocal(hitpos,hitpos_local); // transform into local frame
+ */
+ fCellIdGeoMap[iChId]->GetMatrix()->MasterToLocal(hitpos, hitpos_local);
+ return hitpos_local[0];
+}
+
+double CbmTofEventClusterizer::GetLocalY(CbmTofHit* pHit)
+{
+ Int_t iChId = pHit->GetAddress();
+ Double_t hitpos[3];
+ Double_t hitpos_local[3] = {3 * 0.};
+ hitpos[0] = pHit->GetX();
+ hitpos[1] = pHit->GetY();
+ hitpos[2] = pHit->GetZ();
+ fChannelInfo = fDigiPar->GetCell(iChId);
+ Int_t iCh = CbmTofAddress::GetChannelId(iChId);
+ if (NULL == fChannelInfo) {
+ LOG(fatal) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
+ }
+ //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
+ /*
+ gGeoManager->MasterToLocal(hitpos,hitpos_local); // transform into local frame
+ */
+ fCellIdGeoMap[iChId]->GetMatrix()->MasterToLocal(hitpos, hitpos_local);
+ return hitpos_local[1];
+}
+
+void CbmTofEventClusterizer::MasterToLocal(const int iChId, const double* global, double* local)
+{
+ fCellIdGeoMap[iChId]->GetMatrix()->MasterToLocal(global, local);
+}
+
+static int iLogCal = 0;
+Bool_t CbmTofEventClusterizer::CalibHits()
+{
+ // correct Y position dependence of arrival time
+ for (Int_t iHitInd = 0; iHitInd < fTofHitsColl->GetEntriesFast(); iHitInd++) {
+ CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(iHitInd);
+ if (NULL == pHit) continue;
+
+ Int_t iDetId = (pHit->GetAddress() & DetMask);
+ Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
+ if (iSmType == 5 || iSmType == 8) continue; //Diamond or Pads
+ Int_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ Int_t iSm = CbmTofAddress::GetSmId(iDetId);
+ Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
+ int iRpcInd = iSm * iNbRpc + iRpc;
+
+ double dY = GetLocalY(pHit);
+ double fdYRange = fvCPTOffYRange[iSmType][iRpcInd];
+ LOG(debug) << "TSR " << iSmType << iSm << iRpc << ": Y " << dY << ", " << fdYRange;
+ if (fdYRange == 0) continue;
+
+ double fdBinDy = fvCPTOffYBinWidth[iSmType][iRpcInd];
+ double dBin = (dY + fdYRange) / fdBinDy;
+ int iBin = floor(dBin);
+ int iBin1 = ceil(dBin);
+ double dRes = dBin - iBin;
+ //interpolate
+ double dTcor = fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][iBin];
+ dTcor += dRes * (fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][iBin1] - fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][iBin]);
+ if (iLogCal < 10) {
+ iLogCal++;
+ LOG(info) << "TSR " << iSmType << iSm << iRpc << ": Y " << dY << ", B " << dBin << ", " << iBin << ", " << iBin1
+ << ", " << dRes << ", TcorY " << dTcor;
+ }
+
+ pHit->SetTime(pHit->GetTime() - dTcor);
+ }
+ return kTRUE;
+}
diff --git a/reco/detectors/tof/CbmTofEventClusterizer.h b/reco/detectors/tof/CbmTofEventClusterizer.h
index b53f4d8aa2798e20db53d00f0c86c5147a9e07c1..ffa54badcf084bc854fa557bde9435a611f4925f 100644
--- a/reco/detectors/tof/CbmTofEventClusterizer.h
+++ b/reco/detectors/tof/CbmTofEventClusterizer.h
@@ -30,10 +30,15 @@ class CbmTofCell;
class CbmTofFindTracks;
class CbmDigiManager;
class CbmTofCalibrator;
+class CbmTsEventHeader;
+class CbmTimeSlice;
+
#include "CbmMatch.h"
#include "CbmTofAddress.h" // in cbmdata/tof
#include "CbmTofDigi.h"
+#include "TGeoManager.h"
+class TimesliceMetaData;
class TTofCalibData;
class TTrbHeader;
class FairEventHeader;
@@ -140,6 +145,8 @@ public:
inline void PosYMaxScal(Double_t val) { fPosYMaxScal = val; }
inline void SetTotMax(Double_t val) { fTotMax = val; }
inline void SetTotMin(Double_t val) { fTotMin = val; }
+ inline double GetTotMax() { return fTotMax; }
+ inline double GetTotMin() { return fTotMin; }
inline void SetTotMean(Double_t val) { fTotMean = val; }
inline void SetDelTofMax(Double_t val) { fdDelTofMax = val; }
inline void SetTotPreRange(Double_t val) { fTotPreRange = val; }
@@ -152,6 +159,10 @@ public:
inline void SetSel2Sm(Int_t ival) { fSel2Sm = ival; }
inline void SetSel2Rpc(Int_t ival) { fSel2Rpc = ival; }
inline void SetSel2MulMax(Int_t Id) { fSel2MulMax = Id; }
+ inline void SetCluSizeMin(Int_t iVal) { fiCluSizeMin = iVal; }
+ inline int GetCluSizeMin() { return fiCluSizeMin; }
+ inline void SetNbCalHitMin(Int_t iVal) { fNbCalHitMin = iVal; }
+ inline int GetNbCalHitMin() { return fNbCalHitMin; }
inline void SetOutHstFileName(TString OutHstFileName) { fOutHstFileName = OutHstFileName; }
inline void SetCalParFileName(TString CalParFileName) { fCalParFileName = CalParFileName; }
@@ -164,8 +175,12 @@ public:
inline int GetNbHits() { return fiNbHits; }
CbmTofHit* GetHitPointer(int iHit);
CbmMatch* GetMatchPointer(int iHit);
+ CbmMatch* GetMatchIndexPointer(int idx);
inline double GetTotMean() { return fTotMean; }
inline int GetBeamAddr() { return fiBeamRefAddr; }
+ double GetLocalX(CbmTofHit* pHit);
+ double GetLocalY(CbmTofHit* pHit);
+ void MasterToLocal(const Int_t CellId, const Double_t* master, Double_t* local);
//static Double_t f1_xboxe(double *x, double *par); // Fit function
virtual void fit_ybox(const char* hname); // Fit
@@ -184,7 +199,18 @@ public:
void SetWriteDigisInOut(Bool_t bDigis) { fbWriteDigisInOut = bDigis; }
void SetWriteHitsInOut(Bool_t bHits) { fbWriteHitsInOut = bHits; }
void SetAlternativeBranchNames(Bool_t bNames) { fbAlternativeBranchNames = bNames; }
- void SetDeadStrips(Int_t iDet, Int_t ival);
+ void SetDeadStrips(Int_t iDet, UInt_t ival);
+
+ inline void SetEdgeThr(double val) { fdEdgeThr = val; }
+ inline void SetEdgeLen(double val) { fdEdgeLen = val; }
+ inline void SetEdgeTbias(double val) { fdEdgeTbias = val; }
+ inline void SetEdgeFrange(double val) { fdEdgeFrange = val; }
+ inline double GetEdgeThr() { return fdEdgeThr; }
+ inline double GetEdgeLen() { return fdEdgeLen; }
+ inline double GetEdgeTbias() { return fdEdgeTbias; }
+ inline double GetEdgeFrange() { return fdEdgeFrange; }
+ inline void SetModifySigvel(Double_t val) { fdModifySigvel = val; }
+ inline double GetModifySigvel() { return fdModifySigvel; }
protected:
private:
@@ -238,6 +264,7 @@ private:
Bool_t BuildHits();
Bool_t CalibRawDigis();
Bool_t InspectRawDigis();
+ Bool_t CalibHits();
// ToF geometry variables
CbmTofGeoHandler* fGeoHandler;
@@ -248,6 +275,8 @@ private:
TTrbHeader* fTrbHeader;
FairEventHeader* fEvtHeader;
+ const CbmTsEventHeader* fTsHeader;
+ const CbmTimeSlice* fTimeSlice;
// Input variables
const std::vector<CbmMatch>* fTofDigiPointMatches = nullptr; // TOF MC point matches
@@ -293,6 +322,7 @@ private:
// Histograms
TH1* fhClustBuildTime;
+ TH2* fhClustHitsDigi;
TH1* fhHitsPerTracks;
TH1* fhPtsPerHit;
TH1* fhTimeResSingHits;
@@ -323,6 +353,10 @@ private:
std::vector<TH2*> fhRpcDigiDTFD; //[nbDet]
std::vector<TH2*> fhRpcDigiDTMul; //[nbDet]
std::vector<TH1*> fhRpcDigiRate; //[nbDet]
+ std::vector<TH2*> fhRpcDigiTotLeft; //[nbDet]
+ std::vector<TH2*> fhRpcDigiTotRight; //[nbDet]
+ std::vector<TH2*> fhRpcDigiTotDiff; //[nbDet]
+ std::vector<TH2*> fhRpcDigiTotMap; //[nbDet]
std::vector<TH1*> fhRpcCluMul; //[nbDet]
std::vector<TH1*> fhRpcCluRate; //[nbDet]
std::vector<TH1*> fhRpcCluRate10s; //[nbDet]
@@ -379,9 +413,12 @@ private:
std::vector<std::vector<std::vector<std::vector<Double_t>>>> fvCPTotOff; //[nSMT][nRpc][nCh][nbSide]
std::vector<std::vector<std::vector<std::vector<std::vector<Double_t>>>>>
fvCPWalk; //[nSMT][nRpc][nCh][nbSide][nbWalkBins]
+ std::vector<std::vector<std::vector<Double_t>>> fvCPTOffY; //[nSMT][nRpc][nBin]
+ std::vector<std::vector<Double_t>> fvCPTOffYBinWidth; //[nSMT][nRpc]
+ std::vector<std::vector<Double_t>> fvCPTOffYRange; //[nSMT][nRpc]
std::vector<std::vector<std::vector<std::vector<std::list<CbmTofHit*>>>>> fvLastHits; //[nSMT[nSm][nRpc][nCh][NHits]
- std::vector<Int_t> fvDeadStrips; //[nbDet]
+ std::vector<UInt_t> fvDeadStrips; //[nbDet]
std::vector<std::vector<Double_t>> fvTimeLastDigi; //[nbDet][nChannel*2]
std::vector<std::vector<Double_t>> fvTimeFirstDigi; //[nbDet][nChannel*2]
std::vector<std::vector<Double_t>> fvMulDigi; //[nbDet][nChannel*2]
@@ -425,8 +462,10 @@ private:
Int_t fSel2Rpc;
Int_t fSel2Addr;
Int_t fSel2MulMax;
-
+ Int_t fiCluSizeMin;
+ Int_t fNbCalHitMin;
std::map<UInt_t, UInt_t> fDetIdIndexMap;
+ std::map<UInt_t, TGeoPhysicalNode*> fCellIdGeoMap;
std::vector<Int_t> fviDetId;
Double_t fPosYMaxScal;
@@ -461,6 +500,11 @@ private:
Double_t fdMaxTimeDist; // Isn't this just a local variable? Why make it global and preset?!?
Double_t fdMaxSpaceDist; // Isn't this just a local variable? Why make it global and preset?!?
+ Double_t fdEdgeThr; // for Calibrator
+ Double_t fdEdgeLen; // for Calibrator
+ Double_t fdEdgeTbias;
+ Double_t fdEdgeFrange;
+ Double_t fdModifySigvel;
Double_t fdEvent;
Double_t fdStartAnalysisTime;
diff --git a/reco/detectors/tof/CbmTofFindTracks.cxx b/reco/detectors/tof/CbmTofFindTracks.cxx
index e3d66dfb33ee08122d47712f382c3e0349fd032b..cb9e7158bc661a3d8c33bbad0253a11e7bb7474e 100644
--- a/reco/detectors/tof/CbmTofFindTracks.cxx
+++ b/reco/detectors/tof/CbmTofFindTracks.cxx
@@ -132,12 +132,22 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fhTrklXY0_0(NULL)
, fhTrklXY0_1(NULL)
, fhTrklXY0_2(NULL)
+ , fhTrklX0_TX(NULL)
+ , fhTrklY0_TX(NULL)
+ , fhTrklX0_TY(NULL)
+ , fhTrklY0_TY(NULL)
, vhPullX()
, vhPullY()
, vhPullZ()
, vhPullT()
, vhPullTB()
, vhTrefRms()
+ , vhFitDT0()
+ , vhFitT0Err()
+ , vhFitTt()
+ , vhFitTtErr()
+ , vhFitDTMean()
+ , vhFitDTMeanErr()
, vhResidualTBWalk()
, vhResidualYWalk()
, vhXY_AllTracks()
@@ -154,7 +164,7 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fhVTXNorm(NULL)
, fhVTX_XY0(NULL)
, fhVTX_DT0_Norm(NULL)
- , fiStationStatus()
+ , fiStationStatus(0)
, fOutHstFileName("")
, fCalParFileName("")
, fCalOutFileName("./tofFindTracks.hst.root")
@@ -174,11 +184,13 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fhTOff_Smt(NULL)
, fhTOff_Smt_Off(NULL)
, fhDeltaTt_Smt(NULL)
+ , fhDeltaTc_Smt(NULL)
, fhTOff_HMul2(NULL)
, fiCorMode(0)
, fiBeamCounter(-1)
, fiStationMaxHMul(1000)
, fTtTarg(30.)
+ , fTtLight(0.)
, fdTOffScal(1.)
, fVTXNorm(0.)
, fVTX_T(0.)
@@ -254,9 +266,9 @@ InitStatus CbmTofFindTracks::Init()
fEventsColl = dynamic_cast<TClonesArray*>(ioman->GetObject("Event"));
if (!fEventsColl) fEventsColl = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
if (!fEventsColl) { LOG(info) << "CbmEvent not found in input file, assume eventwise input"; }
- else {
- fTofHitArray = new TClonesArray("CbmTofHit");
- }
+
+ fTofHitArray = new TClonesArray("CbmTofHit");
+
// Get TOF hit Array
fTofHitArrayIn = (TClonesArray*) ioman->GetObject("TofHit");
@@ -280,14 +292,14 @@ InitStatus CbmTofFindTracks::Init()
if (fEventsColl) {
fTrackArrayOut = new TClonesArray("CbmTofTracklet", 100);
fTofHitArrayOut = new TClonesArray("CbmTofHit", 100);
- ioman->Register("TofTracks", "TOF", fTrackArrayOut, kTRUE);
- ioman->Register("TofCalHit", "TOF", fTofHitArrayOut, kTRUE);
+ ioman->Register("TofTracks", "TOF", fTrackArrayOut, IsOutputBranchPersistent("TofTracks"));
+ ioman->Register("TofCalHit", "TOF", fTofHitArrayOut, IsOutputBranchPersistent("TofCalHit"));
}
else {
- ioman->Register("TofTracks", "TOF", fTrackArray, kTRUE);
+ ioman->Register("TofTracks", "TOF", fTrackArray, IsOutputBranchPersistent("TofTracks"));
cout << "-I- CbmTofFindTracks::Init:TofTrack array registered" << endl;
-
- ioman->Register("TofCalHit", "TOF", fTofHitArray, kFALSE);
+ LOG(info) << "Register TofCalHit at " << fTofHitArray;
+ ioman->Register("TofCalHit", "TOF", fTofHitArray, IsOutputBranchPersistent("TofCalHit"));
// Create and register TofUHit array for unused Hits
ioman->Register("TofUHit", "TOF", fTofUHitArray, kFALSE);
@@ -488,7 +500,7 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
new TH1F(Form("hPullT_Smt_Off"), Form("Tracklet PullT vs RpcInd ; RpcInd ; #DeltaT (ns)"), nSmt, 0, nSmt);
}
// Initialize Parameter
- if (fiCorMode <= 3) { // hidden option, FIXME
+ if (fiCorMode <= -3) { // hidden option, FIXME, disabled
for (Int_t iDet = 0; iDet < nSmt; iDet++) {
std::map<Int_t, Int_t>::iterator it;
//it = fMapRpcIdParInd.find(iDet)
@@ -506,11 +518,13 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
if (dVal == 0) {
if (fiBeamCounter != iUniqueId) dVal = fChannelInfo->GetZ() * fTtTarg; // use calibration target value
fhPullT_Smt_Off->SetBinContent(iDet + 1, dVal);
+ LOG(info) << Form("PrimInit det 0x%08x at %d, GloInd %d, z=%f to Tt %6.4f, %6.4f with TOff %6.2f", iUniqueId,
+ iDet, iMap, fChannelInfo->GetZ(), fTtTarg, fdTOffScal, dVal);
}
else {
if (fdTOffScal != 0.) fhPullT_Smt_Off->SetBinContent(iDet + 1, dVal * fdTOffScal);
}
- LOG(info) << Form("Initialize det 0x%08x at %d, GloInd %d, z=%f to Tt %6.4f, %6.4f with TOff %6.2f", iUniqueId,
+ LOG(info) << Form("ReInit det 0x%08x at %d, GloInd %d, z=%f to Tt %6.4f, %6.4f with TOff %6.2f", iUniqueId,
iDet, iMap, fChannelInfo->GetZ(), fTtTarg, fdTOffScal, dVal);
}
}
@@ -835,10 +849,12 @@ Bool_t CbmTofFindTracks::WriteHistos()
if (dRMS > fSIGX * 3.0) dRMS = fSIGX * 3.;
// limit maximal shift in X, for larger values, change geometry file
- if (dVal < -3.) dVal = -3.;
- if (dVal > 3.) dVal = 3.;
+ const double dMaxShift = 4.;
+ if (dVal < -dMaxShift) dVal = -dMaxShift;
+ if (dVal > dMaxShift) dVal = dMaxShift;
LOG(info) << "Update hPullX_Smt_Off " << ix << ": " << fhPullX_Smt_Off->GetBinContent(ix + 1) << " + "
- << htmp1D->GetBinContent(ix + 1) << " -> " << dVal << ", Width " << dRMS;
+ << htmp1D->GetBinContent(ix + 1) << ", " << fRes->Parameter(1) << " -> " << dVal << ", Width "
+ << dRMS;
if (fRpcAddr[ix] != fiBeamCounter) // don't correct beam counter position
fhPullX_Smt_Off->SetBinContent(ix + 1, dVal);
fhPullX_Smt_Width->SetBinContent(ix + 1, dRMS);
@@ -1145,7 +1161,7 @@ void CbmTofFindTracks::Exec(Option_t* opt)
if (fTofHitArray) fTofHitArray->Clear("C");
Int_t iNbHits = 0;
fTofHitIndexArray.resize(tEvent->GetNofData(ECbmDataType::kTofHit));
- for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
+ for (Int_t iHit = 0; iHit < (int) tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHitArrayIn->At(iHitIndex));
fTofHitIndexArray[iNbHits] = iHitIndex;
@@ -1184,6 +1200,7 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
if (NULL != fTofUHitArray) fTofUHitArray->Clear("C");
if (NULL != fTrackArray) fTrackArray->Delete(); // reset
//if (NULL != fTrackArray) fTrackArray->Clear(); // reset
+ CbmTofHit* pBeamHit = NULL;
// recalibrate hits and count trackable hits
for (Int_t iHit = 0; iHit < fTofHitArray->GetEntriesFast(); iHit++) {
@@ -1220,6 +1237,11 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
pHit->SetPosition(hitPos);
TVector3 hitPosErr0(1., 1., 1.); // including positioning uncertainty
pHit->SetPositionError(hitPosErr0); //
+ if (NULL != pBeamHit) {
+ if (pHit->GetTime() < pBeamHit->GetTime()) pBeamHit = pHit;
+ }
+ else
+ pBeamHit = pHit;
}
}
@@ -1260,16 +1282,24 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
LOG(fatal) << "Invalid NTofStations " << fNTofStations << ", " << fMapStationRpcId.size();
}
- LOG(debug) << Form("Exec found Hit %02d, addr 0x%08x, TSR %d%d%d, sta %02d, %02d, HM "
+ LOG(debug) << Form("ExecFind Hit %02d, addr 0x%08x, TSR %d%d%d, sta %02d, %02d, HM "
"%02d, X %6.2f(%3.2f) Y "
- "%6.2f(%3.2f) Z %6.2f(%3.2f) T %6.2f(%3.2f)",
+ "%6.2f(%3.2f) Z %6.2f(%3.2f) T %12.2f(%3.2f)",
iHit, pHit->GetAddress(), CbmTofAddress::GetSmType(iDetId), CbmTofAddress::GetSmId(iDetId),
CbmTofAddress::GetRpcId(iDetId), GetStationOfAddr(iDetId), fDigiBdfPar->GetTrackingStation(pHit),
fStationHMul[GetStationOfAddr(iDetId)], pHit->GetX(), pHit->GetDx(), pHit->GetY(), pHit->GetDy(),
pHit->GetZ(), pHit->GetDz(), pHit->GetTime(), pHit->GetTimeError());
+
MarkStationFired(iSt);
}
+ /* call moved to Clusterizer
+ if (fiCalOpt > 99) {
+ //LOG(info) << "Call TofCalibrator with pBeam "<<pBeamHit;
+ if (NULL != pBeamHit) fTofCalibrator->FillHitCalHist(pBeamHit, fiCalOpt, tEvent, fTofHitArrayIn);
+ }
+ */
+
LOG(debug) << Form("CbmTofFindTracks::Exec NStationsFired %d >= %d Min ?, NbStations %d", GetNStationsFired(),
GetMinNofHits(), fDigiBdfPar->GetNbTrackingStations());
@@ -1308,10 +1338,11 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
void CbmTofFindTracks::Finish()
{
if (fiEvent < 1000) return; // preserve calibration histos in event display
+ FairLogger::GetLogger()->SetLogScreenLevel("info");
if (fiCalOpt > 0) fTofCalibrator->UpdateCalHist(fiCalOpt);
if (fbDoHistos) WriteHistos();
- LOG(info) << Form(" CbmTofFindTracks::Finished ");
+ LOG(info) << Form(" CbmTofFindTracks::Finished for opt %d", fiCalOpt);
}
// -------------------------------------------------------------------------
@@ -1393,10 +1424,19 @@ void CbmTofFindTracks::CreateHistograms()
Double_t X0MAX = 40.;
fhTrklXY0_0 = new TH2F(Form("hTrklXY0_0"), Form("Tracklet XY at z=0 for hmulmax ; x (cm); y (cm)"), 100, -X0MAX,
X0MAX, 100, -X0MAX, X0MAX);
- fhTrklXY0_1 = new TH2F(Form("hTrklXY0_1"), Form("Tracklet XY at z=0 for hmulmax-1 ; x (cm); y (cm)"), 100,
- -X0MAX * 2., X0MAX * 2., 100, -X0MAX * 2., X0MAX * 2.);
- fhTrklXY0_2 = new TH2F(Form("hTrklXY0_2"), Form("Tracklet XY at z=0 for hmulmax-2 ; x (cm); y (cm)"), 100,
- -X0MAX * 3., X0MAX * 3., 100, -X0MAX * 3., X0MAX * 3.);
+ fhTrklXY0_1 = new TH2F(Form("hTrklXY0_1"), Form("Tracklet XY at z=0 for hmulmax-1 ; x (cm); y (cm)"), 100, -X0MAX,
+ X0MAX, 100, -X0MAX, X0MAX);
+ fhTrklXY0_2 = new TH2F(Form("hTrklXY0_2"), Form("Tracklet XY at z=0 for hmulmax-2 ; x (cm); y (cm)"), 100, -X0MAX,
+ X0MAX, 100, -X0MAX, X0MAX);
+
+ fhTrklX0_TX = new TH2F(Form("hTrklX0_TX"), Form("Tracklet X0 vs TX at z=0 for hmulmax ; tx (); x0 (cm)"), 100, 0.,
+ 0.4, 100, -X0MAX, X0MAX);
+ fhTrklY0_TX = new TH2F(Form("hTrklY0_TX"), Form("Tracklet Y0 vs TX at z=0 for hmulmax ; tx (); y0 (cm)"), 100, 0.,
+ 0.4, 100, -X0MAX, X0MAX);
+ fhTrklX0_TY = new TH2F(Form("hTrklX0_TY"), Form("Tracklet X0 vs TY at z=0 for hmulmax ; ty (); x0 (cm)"), 100, 0.,
+ 0.4, 100, -X0MAX, X0MAX);
+ fhTrklY0_TY = new TH2F(Form("hTrklY0_TY"), Form("Tracklet Y0 vs TY at z=0 for hmulmax ; ty (); y0 (cm)"), 100, 0.,
+ 0.4, 100, -X0MAX, X0MAX);
Double_t DT0MAX = 5.;
if (fT0MAX == 0) fT0MAX = DT0MAX;
@@ -1421,12 +1461,21 @@ void CbmTofFindTracks::CreateHistograms()
fhDeltaTt_Smt = new TH2F(Form("hDeltaTt_Smt"), Form("Tracklet DeltaTt; RpcInd ; #DeltaTt (ns/cm)"), nSmt, 0, nSmt,
100, -DTTMAX, DTTMAX);
+ fhDeltaTc_Smt =
+ new TH2F(Form("hDeltaTc_Smt"), Form("Tracklet DeltaTc; RpcInd ; #DeltaTc (ns)"), nSmt, 0, nSmt, 160, -2., 2.);
+
vhPullX.resize(fNTofStations);
vhPullY.resize(fNTofStations);
vhPullZ.resize(fNTofStations);
vhPullT.resize(fNTofStations);
vhPullTB.resize(fNTofStations);
vhTrefRms.resize(fNTofStations);
+ vhFitDT0.resize(fNTofStations);
+ vhFitT0Err.resize(fNTofStations);
+ vhFitTt.resize(fNTofStations);
+ vhFitTtErr.resize(fNTofStations);
+ vhFitDTMean.resize(fNTofStations);
+ vhFitDTMeanErr.resize(fNTofStations);
vhResidualTBWalk.resize(fNTofStations);
vhResidualYWalk.resize(fNTofStations);
vhXY_AllTracks.resize(fNTofStations);
@@ -1453,6 +1502,19 @@ void CbmTofFindTracks::CreateHistograms()
vhPullTB[iSt] = new TH1F(Form("hPullTB_Station_%d", iSt), Form("hResiTB_Station_%d; #DeltaT (ns)", iSt), 59,
-1.25 * fT0MAX, 1.25 * fT0MAX);
vhTrefRms[iSt] = new TH1F(Form("hTrefRms_Station_%d", iSt), Form("hTrefRms_Station_%d; RMS (ns)", iSt), 40, 0, 2.);
+ vhFitDT0[iSt] =
+ new TH1F(Form("hFitDT0_Station_%d", iSt), Form("hFitDT0_Station_%d; #Deltat (ns)", iSt), 40, -5, 5.);
+ vhFitT0Err[iSt] =
+ new TH1F(Form("hFitT0Err_Station_%d", iSt), Form("hFitT0Err_Station_%d; rel. error ()", iSt), 40, 0, 0.1);
+ vhFitTt[iSt] =
+ new TH1F(Form("hFitTt_Station_%d", iSt), Form("hFitTt_Station_%d; inverse velocity (ns/cm)", iSt), 40, 0, 0.1);
+ vhFitTtErr[iSt] = new TH1F(Form("hFitTtErr_Station_%d", iSt),
+ Form("hFitTtErr_Station_%d; inverse velocity error()", iSt), 40, 0, 0.4);
+ vhFitDTMean[iSt] = new TH1F(Form("hFitDTMean_Station_%d", iSt),
+ Form("hFitDTMean_Station_%d; extrapolation distance (ns)", iSt), 40, -5, 5.);
+ vhFitDTMeanErr[iSt] = new TH1F(Form("hFitDTMeanErr_Station_%d", iSt),
+ Form("hFitDTMeanErr_Station_%d; extrapolation error (ns)", iSt), 40, 0, 0.4);
+
const Double_t TOTmax = 50.;
vhResidualTBWalk[iSt] =
new TH2F(Form("hResidualTBWalk_Station_%d", iSt), Form("hResidualTBWalk_Station_%d; #DeltaT (ns)", iSt), TOTmax,
@@ -1464,7 +1526,7 @@ void CbmTofFindTracks::CreateHistograms()
Int_t Nbins = 32.;
CbmTofCell* fChannelInfo = fDigiPar->GetCell(fMapStationRpcId[iSt]);
if (NULL == fChannelInfo) {
- LOG(fatal) << "Geometry for station " << iSt << Form(", Rpc %0x08x ", fMapStationRpcId[iSt]) << " not defined ";
+ LOG(fatal) << "Geometry for station " << iSt << Form(", Rpc 0x%08x ", fMapStationRpcId[iSt]) << " not defined ";
return;
}
Int_t NbinsX = 2 * XSIZ / fChannelInfo->GetSizex();
@@ -1485,7 +1547,7 @@ void CbmTofFindTracks::CreateHistograms()
vhXY_DT[iSt] = new TH3F(Form("hXY_DT_%d", iSt), Form("hXY_DT_%d; x(cm); y (cm); #DeltaT (ns)", iSt), NbinsX, -XSIZ,
XSIZ, Nbins, -XSIZ, XSIZ, Nbins, -0.5, 0.5);
vhXY_TOT[iSt] = new TH3F(Form("hXY_TOT_%d", iSt), Form("hXY_TOT_%d; x(cm); y (cm); TOT (a.u.)", iSt), NbinsX,
- -XSIZ, XSIZ, Nbins, -XSIZ, XSIZ, Nbins, 0., 10.);
+ -XSIZ, XSIZ, Nbins, -XSIZ, XSIZ, Nbins, 0., 20.);
vhXY_CSZ[iSt] = new TH3F(Form("hXY_CSZ_%d", iSt), Form("hXY_CSZ_%d; x(cm); y (cm); CSZ ()", iSt), NbinsX, -XSIZ,
XSIZ, Nbins, -XSIZ, XSIZ, 6, 1., 7.);
vhUDXDY_DT[iSt] = new TH3F(Form("hUDXDY_DT_%d", iSt),
@@ -1618,7 +1680,9 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
iTMul++;
fhTrklChi2->Fill(pTrk->GetNofHits(), pTrk->GetChiSq());
- if (fiCalOpt > 0)
+ fhTrklChi2->Fill(pTrk->GetNofHits(), pTrk->GetChiSq());
+
+ if (fiCalOpt > 0) // && fiCalOpt < 10 )
if (pTrk->GetTt() > fdTtMin) fTofCalibrator->FillCalHist(pTrk, fiCalOpt, tEvent);
CbmTofTrackletParam* tPar = pTrk->GetTrackParameter();
@@ -1648,6 +1712,10 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
switch (GetNReqStations() - pTrk->GetNofHits()) {
case 0: // max hit number
fhTrklXY0_0->Fill(pTrk->GetFitX(0.), pTrk->GetFitY(0.));
+ fhTrklX0_TX->Fill(tPar->GetTx(), pTrk->GetFitX(0.));
+ fhTrklY0_TX->Fill(tPar->GetTx(), pTrk->GetFitY(0.));
+ fhTrklX0_TY->Fill(tPar->GetTy(), pTrk->GetFitX(0.));
+ fhTrklY0_TY->Fill(tPar->GetTy(), pTrk->GetFitY(0.));
break;
case 1: fhTrklXY0_1->Fill(pTrk->GetFitX(0.), pTrk->GetFitY(0.)); break;
case 2: fhTrklXY0_2->Fill(pTrk->GetFitX(0.), pTrk->GetFitY(0.)); break;
@@ -1676,9 +1744,9 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
Double_t dDX = pHit->GetX() - pTrk->GetFitX(pHit->GetZ()); // - tPar->GetX() - tPar->GetTx()*dDZ;
Double_t dDY = pHit->GetY() - pTrk->GetFitY(pHit->GetZ()); // - tPar->GetTy()*dDZ;
Double_t dDT = pHit->GetTime() - pTrk->GetFitT(pHit->GetZ()); // pTrk->GetTdif(fMapStationRpcId[iSt]);
- Double_t dTexp = fTrackletTools->GetTexpected(pTrk, fMapStationRpcId[iSt], pHit);
- Double_t dDTB = pHit->GetTime() - dTexp;
+ Double_t dTexp = fTrackletTools->GetTexpected(pTrk, fMapStationRpcId[iSt], pHit, fTtLight);
// fTrackletTools->GetTdif(pTrk, fMapStationRpcId[iSt],pHit); // ignore pHit in calc of reference
+ Double_t dDTB = pHit->GetTime() - dTexp;
Double_t dTexpErr = fTrackletTools->GetTexpectedError(pTrk, fMapStationRpcId[iSt], pHit, dTexp);
Double_t dTOT = pHit->GetCh() / 10.; // misuse of channel field
@@ -1694,12 +1762,33 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
vhPullT[iSt]->Fill(dDT);
vhPullTB[iSt]->Fill(dDTB);
vhTrefRms[iSt]->Fill(dTexpErr);
+ vhFitDT0[iSt]->Fill(pHit->GetTime() - pTrk->GetT0());
+ vhFitT0Err[iSt]->Fill(pTrk->GetT0Err());
+ vhFitTt[iSt]->Fill(pTrk->GetTt());
+ vhFitTtErr[iSt]->Fill(pTrk->GetTtErr() / pTrk->GetTt());
+ vhFitDTMean[iSt]->Fill(fTrackletTools->GetDTMean(pTrk, pHit));
+ vhFitDTMeanErr[iSt]->Fill(fTrackletTools->GetDTMeanError(pTrk, pHit));
vhResidualTBWalk[iSt]->Fill(dTOT, dDTB);
vhResidualYWalk[iSt]->Fill(dTOT, dDY);
fhPullT_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDT);
fhPullX_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDX);
fhPullY_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDY);
+ bool bCalModule = kTRUE;
+ if (bCalModule) { /* fill all counters of module*/
+ int iModType = CbmTofAddress::GetSmType(fMapStationRpcId[iSt]);
+ if (iModType < 3) {
+ int iSm = CbmTofAddress::GetSmId(fMapStationRpcId[iSt]);
+ int iRpc = CbmTofAddress::GetRpcId(fMapStationRpcId[iSt]);
+ for (int i = 0; i < fDigiBdfPar->GetNbDet(); i++) {
+ if (i != iRpc) {
+ int iRpcId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iModType);
+ fhPullX_Smt->Fill((Double_t) fMapRpcIdParInd[iRpcId], dDX);
+ fhPullY_Smt->Fill((Double_t) fMapRpcIdParInd[iRpcId], dDY);
+ }
+ }
+ }
+ }
/*
fhPullT_Smt->Fill((Double_t)fMapRpcIdParInd[fMapStationRpcId[iSt]], fTrackletTools->GetTdif(pTrk,fMapStationRpcId[iSt], pHit) );
fhPullX_Smt->Fill((Double_t)fMapRpcIdParInd[fMapStationRpcId[iSt]], fTrackletTools->GetXdif(pTrk,fMapStationRpcId[iSt], pHit) );
@@ -1709,6 +1798,10 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
Double_t dDeltaTt = dTt - fTtTarg;
fhDeltaTt_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDeltaTt);
+ if (fiBeamCounter != -1) {
+ Double_t dDeltaTc = pHit->GetTime() - pRefHit->GetTime() - pHit->GetR() / 29.979;
+ fhDeltaTc_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDeltaTc);
+ }
//XXX use BRef as Referenz!!!
if (pRefHit != NULL) {
if (pHit != pRefHit) {
@@ -1776,7 +1869,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
hitpos[0] = pHit->GetX() - (Double_t) fhPullX_Smt_Off->GetBinContent(iRpcInd + 1);
hitpos[1] = pHit->GetY() - (Double_t) fhPullY_Smt_Off->GetBinContent(iRpcInd + 1);
hitpos[2] = pHit->GetZ() - (Double_t) fhPullZ_Smt_Off->GetBinContent(iRpcInd + 1);
- /* TGeoNode* cNode= gGeoManager->GetCurrentNode();*/
+ /* TGeoNode* cNode= */ gGeoManager->GetCurrentNode();
gGeoManager->MasterToLocal(hitpos, hitpos_local);
}
vhXY_AllStations[iSt]->Fill(hitpos_local[0], hitpos_local[1]);
@@ -2286,3 +2379,17 @@ int CbmTofFindTracks::GetStationStatus(int iStation)
else
return 0;
}
+
+int CbmTofFindTracks::GetNbHits() { return fTofHitArray->GetEntriesFast(); } // for internal use by Calibrator
+
+CbmTofHit* CbmTofFindTracks::GetHitPointer(int iHit) { return (CbmTofHit*) fTofHitArray->At(iHit); }
+
+int CbmTofFindTracks::GetHitIndex(int iHit)
+{
+ if ((int) fTofHitIndexArray.size() > iHit) {
+ //LOG(debug) << "Hit " << iHit << ", index " << fTofHitIndexArray[iHit];
+ return fTofHitIndexArray[iHit];
+ }
+ else
+ return iHit;
+}
diff --git a/reco/detectors/tof/CbmTofFindTracks.h b/reco/detectors/tof/CbmTofFindTracks.h
index 7a6d799fbc0e005bbd227158f34fd9fbd0147653..6633f84380f37e50944ae065fa76de1080713510 100644
--- a/reco/detectors/tof/CbmTofFindTracks.h
+++ b/reco/detectors/tof/CbmTofFindTracks.h
@@ -56,7 +56,7 @@ public:
CbmTofFindTracks();
- /** Standard constructor
+ /** Standard constructor
**
*@param name Name of class
*@param title Task title
@@ -132,6 +132,7 @@ public:
inline Int_t GetBeamCounter() const { return fiBeamCounter; }
inline Int_t GetEventNumber() const { return fiEvent; }
inline Double_t GetTtTarg() const { return fTtTarg; }
+ inline Double_t GetTtLight() const { return fTtLight; }
inline Double_t GetTOffScal() const { return fdTOffScal; }
inline Double_t GetSigT() const { return fSIGT; }
@@ -163,6 +164,8 @@ public:
inline void SetCalParFileName(TString CalParFileName) { fCalParFileName = CalParFileName; }
inline void SetCalOutFileName(TString CalOutFileName) { fCalOutFileName = CalOutFileName; }
inline void SetTtTarg(Double_t val) { fTtTarg = val; }
+ inline void SetTtLight(Double_t val) { fTtLight = val; }
+
inline void SetTOffScal(Double_t val) { fdTOffScal = val; }
inline void SetT0MAX(Double_t val) { fT0MAX = val; }
@@ -190,6 +193,12 @@ public:
else
return fTofHitIndexArray[iHit];
}
+
+ int GetNbHits();
+ CbmTofHit* GetHitPointer(int iHit);
+ int GetHitIndex(int iHit);
+
+
bool EvalDoublets(int iI0, int iI1, int iI2, double* dTshift);
private:
@@ -261,6 +270,10 @@ private:
TH2* fhTrklXY0_0;
TH2* fhTrklXY0_1;
TH2* fhTrklXY0_2;
+ TH2* fhTrklX0_TX;
+ TH2* fhTrklY0_TX;
+ TH2* fhTrklX0_TY;
+ TH2* fhTrklY0_TY;
std::vector<TH1*> vhPullX;
std::vector<TH1*> vhPullY;
@@ -268,6 +281,12 @@ private:
std::vector<TH1*> vhPullT;
std::vector<TH1*> vhPullTB;
std::vector<TH1*> vhTrefRms;
+ std::vector<TH1*> vhFitDT0;
+ std::vector<TH1*> vhFitT0Err;
+ std::vector<TH1*> vhFitTt;
+ std::vector<TH1*> vhFitTtErr;
+ std::vector<TH1*> vhFitDTMean;
+ std::vector<TH1*> vhFitDTMeanErr;
std::vector<TH2*> vhResidualTBWalk;
std::vector<TH2*> vhResidualYWalk;
std::vector<TH2*> vhXY_AllTracks; // for monitoring
@@ -311,11 +330,13 @@ private:
TH2* fhTOff_Smt; // Time calibration histo
TH1* fhTOff_Smt_Off; // Time calibration histo
TH2* fhDeltaTt_Smt; // Time calibration histo
+ TH2* fhDeltaTc_Smt; // Time calibration histo
TH2* fhTOff_HMul2; // Time calibration histo
Int_t fiCorMode;
Int_t fiBeamCounter;
Int_t fiStationMaxHMul;
Double_t fTtTarg; // expected average slope in ns/cm
+ Double_t fTtLight; // slope of Light in ns/cm
Double_t fdTOffScal; // modifier to tune average velocity
Double_t fVTXNorm; // Number of Hits contributing to Vertex determination
Double_t fVTX_T; // measured event wise t0
diff --git a/reco/detectors/tof/CbmTofTrackletTools.cxx b/reco/detectors/tof/CbmTofTrackletTools.cxx
index 43978e982214eeabb9a22bdae821b4d5419410b2..c6aab26d3336823a10f4ab9959f0f9ef72c99c7c 100644
--- a/reco/detectors/tof/CbmTofTrackletTools.cxx
+++ b/reco/detectors/tof/CbmTofTrackletTools.cxx
@@ -69,7 +69,7 @@ Double_t CbmTofTrackletTools::FitTt(CbmTofTracklet* pTrk, Int_t iDetId)
// Call with iDetId=-1 for full Tracklet fit
Int_t nValidHits = 0;
Int_t iHit1 = 0;
- Double_t dDist1 = 100.;
+ Double_t dDist1 = 0.;
//LOG(info) << "FitTt: " << pTrk->GetNofHits() << " hits, exclude " << Form("0x%08x",iDetId);
Double_t aR[pTrk->GetNofHits()];
Double_t at[pTrk->GetNofHits()];
@@ -79,7 +79,7 @@ Double_t CbmTofTrackletTools::FitTt(CbmTofTracklet* pTrk, Int_t iDetId)
if (nValidHits == 0) {
iHit1 = iHit;
//LOG(info) << "Init Dist1 with " << iHit1;
- Double_t* dPar = Line3DFit(pTrk, iDetId); // spatial fit without iDetId
+ Double_t* dPar = Line3DFit(pTrk, iDetId); // spatial fit without iDetId, dPar[1] = slope dx/dz, dPar[3]=dy/dz
dDist1 = pTrk->GetTofHitPointer(iHit1)->GetZ() * TMath::Sqrt(1. + dPar[1] * dPar[1] + dPar[3] * dPar[3]);
//LOG(info) << "Dist1 = " << dDist1;
}
@@ -116,6 +116,7 @@ Double_t CbmTofTrackletTools::FitTt(CbmTofTracklet* pTrk, Int_t iDetId)
Double_t dT0 = (RRsum * tsum - Rsum * Rtsum) / det_cov_mat; // Best Guess for time at origin
Double_t dTt = (-Rsum * tsum + esum * Rtsum) / det_cov_mat; // Best guess for inverted velocity
Double_t dT0Err = TMath::Sqrt(RRsum / det_cov_mat); // sqrt of (0,0) in Covariance matrix -> error on fT0
+ dT0Err /= dT0; // relative error
Double_t dTtErr = TMath::Sqrt(esum / det_cov_mat); // sqrt of (1,1) in Covariance Matrix -> error on fTt
Double_t dT0TtCov = -Rsum / det_cov_mat; // (0,1)=(1,0) in Covariance Matrix -> cov(fT0,fTt)
dT0 += yoffset; // Adding yoffset again
@@ -263,13 +264,14 @@ Double_t CbmTofTrackletTools::GetTdif(CbmTofTracklet* pTrk, Int_t iDetId, CbmTof
return dTdif;
}
-Double_t CbmTofTrackletTools::GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit)
+Double_t CbmTofTrackletTools::GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit, double TtLight)
{
Double_t dTref = 0.;
Double_t Nref = 0;
Double_t dTt = 0.;
-
- dTt = FitTt(pTrk, iDetId);
+ if (TtLight == 0.) dTt = FitTt(pTrk, iDetId);
+ else
+ dTt = TtLight; // fix slope to speed of light
for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit)) continue;
@@ -283,6 +285,7 @@ Double_t CbmTofTrackletTools::GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, C
return 1.E20;
}
dTref /= (Double_t) Nref;
+
return dTref;
}
@@ -319,4 +322,21 @@ Double_t CbmTofTrackletTools::GetTexpectedError(CbmTofTracklet* pTrk, Int_t iDet
return dTrms;
}
+Double_t CbmTofTrackletTools::GetDTMean(CbmTofTracklet* pTrk, CbmTofHit* pHit)
+{
+ double dTmean = 0.;
+ int nHits = 0;
+ for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
+ if (pTrk->GetTofHitPointer(iHit) == pHit) continue;
+ dTmean += pTrk->GetTofHitPointer(iHit)->GetTime();
+ nHits++;
+ }
+ dTmean /= nHits;
+ return pHit->GetTime() - dTmean;
+}
+
+Double_t CbmTofTrackletTools::GetDTMeanError(CbmTofTracklet* pTrk, CbmTofHit* pHit)
+{
+ return pTrk->GetTtErr() / pTrk->GetTt() * TMath::Abs(GetDTMean(pTrk, pHit));
+}
ClassImp(CbmTofTrackletTools)
diff --git a/reco/detectors/tof/CbmTofTrackletTools.h b/reco/detectors/tof/CbmTofTrackletTools.h
index 46c35558287d8d70205f4b93d2a1de6a559fd523..3936a2d8ee998bcd0cc4a072f7159e158e944756 100644
--- a/reco/detectors/tof/CbmTofTrackletTools.h
+++ b/reco/detectors/tof/CbmTofTrackletTools.h
@@ -39,8 +39,10 @@ public:
virtual Double_t GetXdif(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit);
virtual Double_t GetYdif(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit);
virtual Double_t GetTdif(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit);
- virtual Double_t GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit);
+ virtual Double_t GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit, double TtLight = 0.);
virtual Double_t GetTexpectedError(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit, Double_t dTexp);
+ virtual Double_t GetDTMean(CbmTofTracklet* pTrk, CbmTofHit* pHit);
+ virtual Double_t GetDTMeanError(CbmTofTracklet* pTrk, CbmTofHit* pHit);
private:
ClassDef(CbmTofTrackletTools, 1);