diff --git a/reco/detectors/tof/CbmTofExtendTracks.cxx b/reco/detectors/tof/CbmTofExtendTracks.cxx
index 496cd25fc1ad0044c85a81d8eadba08ca21219e6..ea1b341b4311289c7a351a46d9664145faf2d3dd 100644
--- a/reco/detectors/tof/CbmTofExtendTracks.cxx
+++ b/reco/detectors/tof/CbmTofExtendTracks.cxx
@@ -63,6 +63,7 @@ using std::endl;
using std::vector;
const Int_t NDefSetup = 100;
+const double dStDZ = 1.;
static LKFMinuit fMinuit;
static Int_t fiTS = 0;
static Double_t dSUT_z = 0.;
@@ -223,22 +224,22 @@ InitStatus CbmTofExtendTracks::Init()
// Get Sts hit Array
fStsHitArrayIn = (TClonesArray*) ioman->GetObject("StsHit");
if (!fStsHitArrayIn) {
- LOG(fatal) << "-W- CbmTofExtendTracks::Init: No StsHit array!";
- return kERROR;
+ LOG(warn) << "-W- CbmTofExtendTracks::Init: No StsHit array!";
+ //return kERROR;
}
// Get Much hit Array
fMuchHitArrayIn = (TClonesArray*) ioman->GetObject("MuchPixelHit");
if (!fMuchHitArrayIn) {
- LOG(fatal) << "-W- CbmTofExtendTracks::Init: No MuchHit array!";
- return kERROR;
+ LOG(warn) << "-W- CbmTofExtendTracks::Init: No MuchHit array!";
+ //return kERROR;
}
// Get Rich hit Array
fRichHitArrayIn = (TClonesArray*) ioman->GetObject("RichHit");
if (!fRichHitArrayIn) {
- LOG(fatal) << "-W- CbmTofExtendTracks::Init: No RichHit array!";
- return kERROR;
+ LOG(warn) << "-W- CbmTofExtendTracks::Init: No RichHit array!";
+ //return kERROR;
}
if (kFALSE == InitParameters()) return kFATAL;
@@ -579,13 +580,13 @@ void CbmTofExtendTracks::Exec(Option_t* opt)
// ExecExtend(opt);
}
else {
- LOG(info) << "ExtExec TS " << fiTS << " with " << fEventsColl->GetEntriesFast() << " evts";
+ LOG(debug) << "ExtExec TS " << fiTS << " with " << fEventsColl->GetEntriesFast() << " evts";
fiTS++;
for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
- LOG(info) << "Process TS event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kT0Hit) << " T0, "
- << tEvent->GetNofData(ECbmDataType::kTofHit) << " TOF, " << tEvent->GetNofData(ECbmDataType::kStsHit)
- << " STS, " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit) << " MUCH hits";
+ LOG(debug) << "Process TS event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kT0Hit) << " T0, "
+ << tEvent->GetNofData(ECbmDataType::kTofHit) << " TOF, " << tEvent->GetNofData(ECbmDataType::kStsHit)
+ << " STS, " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit) << " MUCH hits";
ExecExtend(opt, tEvent);
}
@@ -601,8 +602,8 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
UInt_t iNbStsStations = 0;
UInt_t iNbMuchStations = 0;
UInt_t iNbRichStations = 0;
- const Double_t STATION_TOF_ZWIDTH = 1.;
- const Double_t STATION_STS_ZWIDTH = 0.5;
+ const Double_t STATION_TOF_ZWIDTH = 10.;
+ const Double_t STATION_STS_ZWIDTH = 1.5;
const Double_t STATION_MUCH_ZWIDTH = 1.;
const Double_t STATION_RICH_ZWIDTH = 3.;
fvTofHitIndex.clear();
@@ -610,11 +611,11 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
fvMuchHitIndex.clear();
fvRichHitIndex.clear();
- LOG(info) << "Ev " << fiEvent << " has hits "
- << " Tof: " << tEvent->GetNofData(ECbmDataType::kTofHit)
- << " Sts: " << tEvent->GetNofData(ECbmDataType::kStsHit)
- << " Much: " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit)
- << " Rich: " << tEvent->GetNofData(ECbmDataType::kRichHit);
+ LOG(debug) << "Ev " << fiEvent << " has hits "
+ << " Tof: " << tEvent->GetNofData(ECbmDataType::kTofHit)
+ << " Sts: " << tEvent->GetNofData(ECbmDataType::kStsHit)
+ << " Much: " << tEvent->GetNofData(ECbmDataType::kMuchPixelHit)
+ << " Rich: " << tEvent->GetNofData(ECbmDataType::kRichHit);
// cleanup
/*
@@ -640,7 +641,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
LOG(debug) << Form("Inspect Ev %d, TofHit %d, Ind %d at %6.1f in %lu (%u) stations", fiEvent, iHit, iHitIndex,
tHit->GetZ(), fvTofStationZ.size(), iNbAllStations);
- Int_t iStZ = (Int_t)(tHit->GetZ());
+ Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
itMapStationZ = fMapStationZ.find(iStZ);
UInt_t iSt = 0;
@@ -679,7 +680,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
}
if (fiEvent < NDefSetup) {
if (itMapStationZ == fMapStationZ.end()) {
- LOG(debug) << "Insert new tracking station " << Int_t(ECbmModuleId::kTof) * 100 + iSt << " at z=" << iStZ;
+ LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kTof) * 100 + iSt << " at z=" << iStZ;
fMapStationZ[iStZ] = Int_t(ECbmModuleId::kTof) * 100 + iSt;
itMapStationZ = fMapStationZ.begin();
Int_t iStId = Int_t(ECbmModuleId::kTof) * 100;
@@ -720,7 +721,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kStsHit, iHit));
CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fStsHitArrayIn->At(iHitIndex));
- Int_t iStZ = (Int_t)(tHit->GetZ());
+ Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
itMapStationZ = fMapStationZ.find(iStZ);
UInt_t iSt = 0;
@@ -758,7 +759,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
}
if (fiEvent < NDefSetup) {
if (itMapStationZ == fMapStationZ.end()) {
- LOG(debug) << "Insert new tracking station " << Int_t(ECbmModuleId::kSts) * 100 + iSt << " at z=" << iStZ;
+ LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kSts) * 100 + iSt << " at z=" << iStZ;
fMapStationZ[iStZ] = Int_t(ECbmModuleId::kSts) * 100 + iSt;
itMapStationZ = fMapStationZ.begin();
Int_t iStId = Int_t(ECbmModuleId::kSts) * 100;
@@ -794,7 +795,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kMuchPixelHit, iHit));
CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fMuchHitArrayIn->At(iHitIndex));
- Int_t iStZ = (Int_t)(tHit->GetZ());
+ Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
itMapStationZ = fMapStationZ.find(iStZ);
UInt_t iSt = 0;
@@ -833,7 +834,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
if (fiEvent < NDefSetup) {
if (itMapStationZ == fMapStationZ.end()) {
- LOG(debug) << "Insert new tracking station " << Int_t(ECbmModuleId::kMuch) * 100 + iSt << " at z=" << iStZ;
+ LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kMuch) * 100 + iSt << " at z=" << iStZ;
fMapStationZ[iStZ] = Int_t(ECbmModuleId::kMuch) * 100 + iSt;
itMapStationZ = fMapStationZ.begin();
Int_t iStId = Int_t(ECbmModuleId::kMuch) * 100;
@@ -886,7 +887,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kRichHit, iHit));
CbmPixelHit* tHit = dynamic_cast<CbmPixelHit*>(fRichHitArrayIn->At(iHitIndex));
- Int_t iStZ = (Int_t)(tHit->GetZ());
+ Int_t iStZ = (Int_t)(tHit->GetZ() / dStDZ);
itMapStationZ = fMapStationZ.find(iStZ);
UInt_t iSt = 0;
@@ -925,7 +926,7 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
if (fiEvent < NDefSetup) {
if (itMapStationZ == fMapStationZ.end()) {
- LOG(debug) << "Insert new tracking station " << Int_t(ECbmModuleId::kRich) * 100 + iSt << " at z=" << iStZ;
+ LOG(info) << "Insert new tracking station " << Int_t(ECbmModuleId::kRich) * 100 + iSt << " at z=" << iStZ;
fMapStationZ[iStZ] = Int_t(ECbmModuleId::kRich) * 100 + iSt;
itMapStationZ = fMapStationZ.begin();
Int_t iStId = Int_t(ECbmModuleId::kRich) * 100;
@@ -989,6 +990,10 @@ void CbmTofExtendTracks::ExecExtend(Option_t* /*opt*/, CbmEvent* tEvent)
for (Int_t iHit = 0; iHit < iNHits; iHit++) {
Int_t iHitInd = pTrk->GetTofHitIndex(iHit);
CbmTofHit* pHitIn = dynamic_cast<CbmTofHit*>(fTofHitArrayIn->At(iHitInd));
+ if (NULL == pHitIn) {
+ LOG(warn) << Form("Hit %d not found at index %d ", iHit, iHitInd);
+ continue;
+ }
//CbmTofHit* pHit = new CbmTofHit(*pHitIn); // copy construction
// apply recalibration if necessary
fvTrkCalHits[iTr].push_back(dynamic_cast<CbmPixelHit*>(pHitIn));
@@ -1140,32 +1145,34 @@ void CbmTofExtendTracks::CreateHistograms()
if (fiStationUT > -1) {
// histos for Station Under Test
itMapStationZ = fMapStationZ.begin();
- for (Int_t iSt = 0; iSt < fiStationUT; iSt++)
+ for (Int_t iSt = 0; iSt < fiStationUT; iSt++) {
itMapStationZ++;
- dSUT_z = itMapStationZ->first;
- LOG(info) << "Create SUT histos for station " << fiStationUT << " at distance " << dSUT_z;
- const Double_t dSUT_RefDx = 0.15;
- const Double_t dSUT_RefDy = 0.3;
- const Double_t dNbinX = 100;
- const Double_t dNbinY = 100;
- const Double_t dNbinZ = 50;
-
- Double_t dSUT_dx = dSUT_RefDx * dSUT_z;
- Double_t dSUT_dy = dSUT_RefDy * dSUT_z;
- fhExtSutXY_Found = new TH2F("hExtSutXY_Found", Form("StationUnderTest %d found hits ; X (cm); Y (cm)", fiStationUT),
- dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
- fhExtSutXY_Missed =
- new TH2F("hExtSutXY_Missed", Form("StationUnderTest %d missed hits ; X (cm); Y (cm)", fiStationUT), dNbinX,
- -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
- fhExtSutXY_DX =
- new TH3F("hExtSutXY_DX", Form("StationUnderTest %d #DeltaX ; X (cm); Y (cm); #DeltaX (cm)", fiStationUT), dNbinX,
- -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
- fhExtSutXY_DY =
- new TH3F("hExtSutXY_DY", Form("StationUnderTest %d #DeltaY ; X (cm); Y (cm); #DeltaY (cm)", fiStationUT), dNbinX,
- -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
- fhExtSutXY_DT =
- new TH3F("hExtSutXY_DT", Form("StationUnderTest %d #DeltaT ; X (cm); Y (cm); #DeltaT (ns)", fiStationUT), dNbinX,
- -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -50., 50.);
+ dSUT_z = itMapStationZ->first;
+ LOG(info) << "Create SUT histos for station " << fiStationUT << " at distance " << dSUT_z;
+ const Double_t dSUT_RefDx = 0.15;
+ const Double_t dSUT_RefDy = 0.3;
+ const Double_t dNbinX = 100;
+ const Double_t dNbinY = 100;
+ const Double_t dNbinZ = 50;
+
+ Double_t dSUT_dx = dSUT_RefDx * dSUT_z;
+ Double_t dSUT_dy = dSUT_RefDy * dSUT_z;
+ fhExtSutXY_Found =
+ new TH2F("hExtSutXY_Found", Form("StationUnderTest %d found hits ; X (cm); Y (cm)", fiStationUT), dNbinX,
+ -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
+ fhExtSutXY_Missed =
+ new TH2F("hExtSutXY_Missed", Form("StationUnderTest %d missed hits ; X (cm); Y (cm)", fiStationUT), dNbinX,
+ -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy);
+ fhExtSutXY_DX =
+ new TH3F("hExtSutXY_DX", Form("StationUnderTest %d #DeltaX ; X (cm); Y (cm); #DeltaX (cm)", fiStationUT),
+ dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
+ fhExtSutXY_DY =
+ new TH3F("hExtSutXY_DY", Form("StationUnderTest %d #DeltaY ; X (cm); Y (cm); #DeltaY (cm)", fiStationUT),
+ dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -10., 10.);
+ fhExtSutXY_DT =
+ new TH3F("hExtSutXY_DT", Form("StationUnderTest %d #DeltaT ; X (cm); Y (cm); #DeltaT (ns)", fiStationUT),
+ dNbinX, -dSUT_dx, dSUT_dx, dNbinY, -dSUT_dy, dSUT_dy, dNbinZ, -50., 50.);
+ }
} // StationUT end
gDirectory->cd(oldir->GetPath()); // <= To prevent histos from being sucked in by the param file of the TRootManager!
@@ -1248,7 +1255,7 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
fvTofHitIndex[iSt][iHit], iSt, fvTofHitIndex.size());
assert(tHit);
}
- Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ())] % 100;
+ Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
if (iTrSt < 0 || iTrSt >= (Int_t) fMapStationZ.size()) {
LOG(error) << "Invalid station index " << iTrSt;
continue;
@@ -1279,7 +1286,7 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
fvStsHitIndex[iSt][iHit], iSt, fvStsHitIndex.size());
assert(tHit);
}
- Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ())] % 100;
+ Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
if (TMath::Abs(dDX) > fdTrkCutDX) continue;
@@ -1306,7 +1313,7 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
fvMuchHitIndex[iSt][iHit], iSt, fvMuchHitIndex.size());
assert(tHit);
}
- Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ())] % 100;
+ Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
LOG(info) << Form("Proc ev %d, Trk %d, St %d, iLe%d MtHit: ", fiEvent, iTr, iTrSt, iLev) << tHit->ToString();
@@ -1346,7 +1353,7 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
fvRichHitIndex[iSt][iHit], iSt, fvRichHitIndex.size());
assert(tHit);
}
- Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ())] % 100;
+ Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
Double_t dDX = tHit->GetX() + fvXoff[iTrSt] - tTrk->GetFitX(tHit->GetZ());
if (TMath::Abs(dDX) > fdTrkCutDX) continue;
@@ -1407,12 +1414,12 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
// LOG(info)<<"CheckII with "<<iReqStations<<","<<iHit;
if (iReqStation != 0 && iHit == 0) continue; //skip this track
}
-
+ /*
LOG(info) << "CheckIII accepted tr " << iTr << " with " << fiReqStations;
for (UInt_t iHit = 0; iHit < fvTrkCalHits[iTr].size(); iHit++) {
LOG(info) << " Hit " << iHit << ", station " << fMapStationZ[(Int_t)(fvTrkCalHits[iTr][iHit]->GetZ())] % 100;
}
-
+ */
fhExt_TrkSizChiSq[iLev]->Fill((Double_t) fvTrkCalHits[iTr].size(), fvTrkPar[iTr]->GetChiSq());
fhExt_TrkSizVel[iLev]->Fill((Double_t) fvTrkCalHits[iTr].size(), 1. / (fvTrkPar[iTr]->GetTt()));
@@ -1422,7 +1429,7 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
Bool_t bSUT_Found = kFALSE;
for (UInt_t iHit = 0; iHit < fvTrkCalHits[iTr].size(); iHit++) {
CbmPixelHit* tHit = fvTrkCalHits[iTr][iHit];
- Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ())] % 100;
+ Int_t iTrSt = fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100;
Double_t dDX = tHit->GetX() - GetFitX(tHit->GetZ(), tTrkPar);
Double_t dDY = tHit->GetY() - GetFitY(tHit->GetZ(), tTrkPar);
@@ -1454,9 +1461,9 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
fhTrkPullDY[0]->Fill(iTrSt, dDY);
fhTrkPullDT[0]->Fill(iTrSt, dDT);
fhTrkStationDXDY[0][iTrSt]->Fill(dDX, dDY);
- LOG(info) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iTrSt, 0, dDX,
- dDY)
- << tHit->ToString();
+ LOG(debug) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iTrSt, 0, dDX,
+ dDY)
+ << tHit->ToString();
} // fvTrkCalHits loop end
if (bSUT_Found == kFALSE) { fhExtSutXY_Missed->Fill(GetFitX(dSUT_z, tTrkPar), GetFitY(dSUT_z, tTrkPar)); }
@@ -1468,11 +1475,11 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
/*
LOG(info)<<"Got AllHit for Tr "<<iTr
<<", St "<<iSt
- <<", StZ "<<fMapStationZ[(Int_t)(tHit->GetZ())]%100
+ <<", StZ "<<fMapStationZ[(Int_t)(tHit->GetZ()/dStDZ)]%100
<<", Hit "<< iHit <<", poi "<< tHit;
LOG(info)<<"tHit: "<<tHit->ToString();
*/
- assert((Int_t) iSt == fMapStationZ[(Int_t)(tHit->GetZ())] % 100);
+ assert((Int_t) iSt == fMapStationZ[(Int_t)(tHit->GetZ() / dStDZ)] % 100);
Double_t dDX = tHit->GetX() - GetFitX(tHit->GetZ(), tTrkPar);
if (TMath::Abs(dDX) > fdTrkCutDX) continue;
@@ -1481,9 +1488,9 @@ void CbmTofExtendTracks::FillHistograms(CbmEvent* tEvent)
Double_t dDT = tHit->GetTime() - GetFitT(tHit->GetZ(), tTrkPar);
if (TMath::Abs(dDT) > fdTrkCutDT) continue;
- LOG(info) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iSt, iLev, dDX,
- dDY)
- << tHit->ToString();
+ LOG(debug) << Form("Proc ev %d, Trk %d, St %d, Lev%d, DX %6.3f, dY %6.3f, MtHit: ", fiEvent, iTr, iSt, iLev,
+ dDX, dDY)
+ << tHit->ToString();
// fill tracking station histos
fhTrkStationDX[iLev]->Fill(iSt, dDX);
@@ -1619,7 +1626,7 @@ void CbmTofExtendTracks::TrkAddStation(Int_t iStation)
CbmPixelHit* tHit = fvAllHitPointer[iStation][MatchHit[iCand]];
/*
LOG(info)<<"Add hit "<<MatchHit[iCand]
- <<", station "<<fMapStationZ[(Int_t)(tHit->GetZ())]%100
+ <<", station "<<fMapStationZ[(Int_t)(tHit->GetZ()/dStDZ)]%100
<<" to track "<<MatchTrk[iCand]
<<" with chi2 "<<MatchChi2Min[iCand];
*/
diff --git a/reco/detectors/tof/CbmTofFindTracks.cxx b/reco/detectors/tof/CbmTofFindTracks.cxx
index 9ffac3e5e70fcd3305f3927fcb82374b7da101a5..e3d66dfb33ee08122d47712f382c3e0349fd032b 100644
--- a/reco/detectors/tof/CbmTofFindTracks.cxx
+++ b/reco/detectors/tof/CbmTofFindTracks.cxx
@@ -60,6 +60,7 @@ using std::vector;
//const Int_t DetMask = 0x3FFFFF; // check for consistency with v14a geometry
const Int_t DetMask = 0x1FFFFF; // check for consistency with v21a geometry
+static int iTS = 0;
ClassImp(CbmTofFindTracks);
@@ -103,6 +104,10 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fvXoff()
, fvYoff()
, fvZoff()
+ , fvTsig()
+ , fvXsig()
+ , fvYsig()
+ , fvZsig()
, fhTrklMul(NULL)
, fhTrklChi2(NULL)
, fhAllHitsStation(NULL)
@@ -117,6 +122,7 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fhTrklZ0yHMul(NULL)
, fhTrklTxHMul(NULL)
, fhTrklTyHMul(NULL)
+ , fhTrklTyTx(NULL)
, fhTrklTtHMul(NULL)
, fhTrklVelHMul(NULL)
, fhTrklT0HMul(NULL)
@@ -131,10 +137,12 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, vhPullZ()
, vhPullT()
, vhPullTB()
+ , vhTrefRms()
, vhResidualTBWalk()
, vhResidualYWalk()
, vhXY_AllTracks()
, vhXY_AllStations()
+ , vhXY_AllFitStations()
, vhXY_MissedStation()
, vhXY_DX()
, vhXY_DY()
@@ -146,6 +154,7 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fhVTXNorm(NULL)
, fhVTX_XY0(NULL)
, fhVTX_DT0_Norm(NULL)
+ , fiStationStatus()
, fOutHstFileName("")
, fCalParFileName("")
, fCalOutFileName("./tofFindTracks.hst.root")
@@ -170,6 +179,7 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fiBeamCounter(-1)
, fiStationMaxHMul(1000)
, fTtTarg(30.)
+ , fdTOffScal(1.)
, fVTXNorm(0.)
, fVTX_T(0.)
, fVTX_X(0.)
@@ -189,13 +199,14 @@ CbmTofFindTracks::CbmTofFindTracks(const char* name, const char* /*title*/, CbmT
, fdRefVelMean(0.)
, fdRefDVel(1.E7)
, fdR0Lim(0.)
+ , fdTtMin(0.)
, fStart()
, fStop()
, fdTrackingTime(0.)
, fdBeamMomentumLab(0.)
, fbRemoveSignalPropagationTime(kFALSE)
, fiBeamMaxHMul(1000)
- , fiCalOpt(0)
+ , fiCalOpt((int) 0)
{
if (!fInstance) fInstance = this;
}
@@ -300,10 +311,8 @@ InitStatus CbmTofFindTracks::Init()
Int_t iCellId = fDigiPar->GetCellId(iCell);
Int_t iCh = fTofId->GetCell(iCellId);
if (0 == iCh) {
- LOG(info) << Form("Init found RpcInd %d, %lu at Addr 0x%08x, ModType %d, "
- "ModId %d, RpcId %d ",
- iRpc, fRpcAddr.size(), iCellId, fTofId->GetSMType(iCellId), fTofId->GetSModule(iCellId),
- fTofId->GetCounter(iCellId));
+ LOG(info) << Form("Init found at Ind %d, %lu Rpc with Addr 0x%08x, TSR %d%d%d ", iRpc, fRpcAddr.size(), iCellId,
+ fTofId->GetSMType(iCellId), fTofId->GetSModule(iCellId), fTofId->GetCounter(iCellId));
if (fTofId->GetSMType(iCellId) == 5) {
bBeamCounter = kTRUE;
LOG(info) << "Found beam counter in setup! at RpcInd " << iRpc << ", Addr.size " << fRpcAddr.size();
@@ -313,8 +322,10 @@ InitStatus CbmTofFindTracks::Init()
iRpc++;
}
}
- fStationHMul.resize(fNTofStations + 1);
+ LOG(debug) << "Initialize fStationHMul to size " << fNTofStations + 1;
+ fStationHMul.resize(fNTofStations + 1);
+ LOG(debug) << "Initialize fStationHMul to size " << fNTofStations + 1;
LoadCalParameter();
@@ -324,7 +335,7 @@ InitStatus CbmTofFindTracks::Init()
fTofCalibrator = new CbmTofCalibrator();
if (fTofCalibrator->Init() != kSUCCESS) return kFATAL;
if (bBeamCounter) {
- fTofCalibrator->SetBeam(bBeamCounter);
+ if (fiBeamCounter > -1) fTofCalibrator->SetBeam(bBeamCounter);
fTofCalibrator->SetR0Lim(fdR0Lim);
LOG(info) << "Set CbmTofCalibrator::R0Lim to " << fdR0Lim;
}
@@ -340,17 +351,23 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
{
UInt_t NSt = fMapRpcIdParInd.size();
fvToff.resize(NSt);
- for (uint i = 0; i < NSt; i++)
- fvToff[i] = 0.;
fvXoff.resize(NSt);
- for (uint i = 0; i < NSt; i++)
- fvXoff[i] = 0.;
fvYoff.resize(NSt);
- for (uint i = 0; i < NSt; i++)
- fvYoff[i] = 0.;
fvZoff.resize(NSt);
- for (uint i = 0; i < NSt; i++)
+ fvTsig.resize(NSt);
+ fvXsig.resize(NSt);
+ fvYsig.resize(NSt);
+ fvZsig.resize(NSt);
+ for (uint i = 0; i < NSt; i++) {
+ fvToff[i] = 0.;
+ fvXoff[i] = 0.;
+ fvYoff[i] = 0.;
fvZoff[i] = 0.;
+ fvTsig[i] = fSIGT;
+ fvXsig[i] = fSIGX;
+ fvYsig[i] = fSIGY;
+ fvZsig[i] = fSIGZ;
+ }
if (fCalParFileName.IsNull()) return kTRUE;
@@ -372,7 +389,7 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
TH1D* fhtmpX = (TH1D*) gDirectory->FindObjectAny(Form("hPullX_Smt_Off"));
TH1D* fhtmpY = (TH1D*) gDirectory->FindObjectAny(Form("hPullY_Smt_Off"));
TH1D* fhtmpZ = (TH1D*) gDirectory->FindObjectAny(Form("hPullZ_Smt_Off"));
- TH1D* fhtmpW = (TH1D*) gDirectory->FindObjectAny(Form("hPullT_Smt_Width"));
+ TH1D* fhtmpWT = (TH1D*) gDirectory->FindObjectAny(Form("hPullT_Smt_Width"));
TH1D* fhtmpWX = (TH1D*) gDirectory->FindObjectAny(Form("hPullX_Smt_Width"));
TH1D* fhtmpWY = (TH1D*) gDirectory->FindObjectAny(Form("hPullY_Smt_Width"));
TH1D* fhtmpWZ = (TH1D*) gDirectory->FindObjectAny(Form("hPullZ_Smt_Width"));
@@ -406,24 +423,60 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
fvZoff[iSt] = fhPullZ_Smt_Off->GetBinContent(iSt + 1);
}
- if (NULL == fhtmpW) { LOG(info) << Form("CbmTofFindTracks::LoadCalParameter: hPullT_Smt_Width") << " not found. "; }
+ if (NULL == fhtmpWT) { LOG(info) << Form("CbmTofFindTracks::LoadCalParameter: hPullT_Smt_Width") << " not found. "; }
else {
- if (fbUseSigCalib) fhPullT_Smt_Width = (TH1D*) fhtmpW->Clone();
+ if (fbUseSigCalib) {
+ fhPullT_Smt_Width = (TH1D*) fhtmpWT->Clone();
+ for (UInt_t iSt = 0; iSt < NSt; iSt++) {
+ fvTsig[iSt] = fhPullT_Smt_Width->GetBinContent(iSt + 1);
+ if (fvTsig[iSt] == 0) {
+ LOG(warning) << "Invalid Tsig for station " << iSt;
+ fvTsig[iSt] = fSIGT;
+ }
+ }
+ }
}
if (NULL == fhtmpWX) { LOG(info) << Form("CbmTofFindTracks::LoadCalParameter: hPullX_Smt_Width") << " not found. "; }
else {
- if (fbUseSigCalib) fhPullX_Smt_Width = (TH1D*) fhtmpWX->Clone();
+ if (fbUseSigCalib) {
+ fhPullX_Smt_Width = (TH1D*) fhtmpWX->Clone();
+ for (UInt_t iSt = 0; iSt < NSt; iSt++) {
+ fvXsig[iSt] = fhPullX_Smt_Width->GetBinContent(iSt + 1);
+ if (fvXsig[iSt] == 0) {
+ LOG(warning) << "Invalid Xsig for station " << iSt;
+ fvXsig[iSt] = fSIGX;
+ }
+ }
+ }
}
if (NULL == fhtmpWY) { LOG(info) << Form("CbmTofFindTracks::LoadCalParameter: hPullY_Smt_Width") << " not found. "; }
else {
- if (fbUseSigCalib) fhPullY_Smt_Width = (TH1D*) fhtmpWY->Clone();
+ if (fbUseSigCalib) {
+ fhPullY_Smt_Width = (TH1D*) fhtmpWY->Clone();
+ for (UInt_t iSt = 0; iSt < NSt; iSt++) {
+ fvYsig[iSt] = fhPullY_Smt_Width->GetBinContent(iSt + 1);
+ if (fvYsig[iSt] == 0) {
+ LOG(warning) << "Invalid Ysig for station " << iSt;
+ fvYsig[iSt] = fSIGY;
+ }
+ }
+ }
}
if (NULL == fhtmpWZ) { LOG(info) << Form("CbmTofFindTracks::LoadCalParameter: hPullZ_Smt_Width") << " not found. "; }
else {
- if (fbUseSigCalib) fhPullZ_Smt_Width = (TH1D*) fhtmpWZ->Clone();
+ if (fbUseSigCalib) {
+ fhPullZ_Smt_Width = (TH1D*) fhtmpWZ->Clone();
+ for (UInt_t iSt = 0; iSt < NSt; iSt++) {
+ fvZsig[iSt] = fhPullZ_Smt_Width->GetBinContent(iSt + 1);
+ if (fvZsig[iSt] == 0) {
+ LOG(warning) << "Invalid Zsig for station " << iSt;
+ fvZsig[iSt] = fSIGZ;
+ }
+ }
+ }
}
fCalParFile->Close();
@@ -433,28 +486,35 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
if (NULL == fhPullT_Smt_Off) { // provide default TOffset histogram
fhPullT_Smt_Off =
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
- for (Int_t iDet = 0; iDet < nSmt; iDet++) {
- std::map<Int_t, Int_t>::iterator it;
- //it = fMapRpcIdParInd.find(iDet);
- for (it = fMapRpcIdParInd.begin(); it != fMapRpcIdParInd.end(); it++) {
- if (it->second == iDet) break;
- }
- LOG(debug1) << Form(" iDet %d -> iUniqueId ? 0x%08x, 0x%08x ", iDet, it->first, it->second);
- Int_t iUniqueId = it->first;
- CbmTofCell* fChannelInfo = fDigiPar->GetCell(iUniqueId);
- if (NULL != fChannelInfo) {
- Double_t dVal = 0.; // FIXME numeric constant in code, default for cosmic
+ }
+ // Initialize Parameter
+ if (fiCorMode <= 3) { // hidden option, FIXME
+ for (Int_t iDet = 0; iDet < nSmt; iDet++) {
+ std::map<Int_t, Int_t>::iterator it;
+ //it = fMapRpcIdParInd.find(iDet)
+ Int_t iMap = 0;
+ for (it = fMapRpcIdParInd.begin(); it != fMapRpcIdParInd.end(); it++) {
+ iMap++;
+ if (it->second == iDet) break;
+ }
+ LOG(debug1) << Form(" iDet %d -> iUniqueId ? 0x%08x, 0x%08x ", iDet, it->first, it->second);
+ Int_t iUniqueId = it->first;
+ CbmTofCell* fChannelInfo = fDigiPar->GetCell(iUniqueId);
+ if (NULL != fChannelInfo) {
+ Double_t dVal = 0.; // FIXME numeric constant in code, default for cosmic
+ dVal = fhPullT_Smt_Off->GetBinContent(iDet + 1, dVal);
+ if (dVal == 0) {
if (fiBeamCounter != iUniqueId) dVal = fChannelInfo->GetZ() * fTtTarg; // use calibration target value
fhPullT_Smt_Off->SetBinContent(iDet + 1, dVal);
- LOG(info) << Form("Initialize det 0x%08x at %d, z=%f with TOff %6.2f", iUniqueId, iDet + 1,
- fChannelInfo->GetZ(), 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,
+ iDet, iMap, fChannelInfo->GetZ(), fTtTarg, fdTOffScal, dVal);
}
+ }
}
-
if (NULL == fhPullT_Smt_Width) { // provide default TWidth histogram
fhPullT_Smt_Width =
new TH1F(Form("hPullT_Smt_Width"), Form("Tracklet ResiT Width vs RpcInd ; RpcInd ; RMS(T) (ns)"), nSmt, 0, nSmt);
@@ -491,7 +551,7 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
}
}
- if (NULL == fhPullZ_Smt_Off) // provide default TOffset histogram
+ if (NULL == fhPullZ_Smt_Off) // provide default ZOffset histogram
fhPullZ_Smt_Off =
new TH1F(Form("hPullZ_Smt_Off"), Form("Tracklet ResiZ vs RpcInd ; RpcInd ; #DeltaZ (cm)"), nSmt, 0, nSmt);
if (NULL == fhPullZ_Smt_Width) {
@@ -678,41 +738,65 @@ Bool_t CbmTofFindTracks::WriteHistos()
{
TProfile* htmp = fhPullT_Smt->ProfileX();
TH1D* htmp1D = htmp->ProjectionX();
-
if (fhPullT_Smt_Off != NULL) {
Double_t nx = htmp1D->GetNbinsX();
for (Int_t ix = 0; ix < nx; ix++) {
- TH1D* hpy = fhPullT_Smt->ProjectionY("_py", ix + 1, ix + 1);
- if (hpy->GetEntries() > 100.) {
+ Double_t dVal = fhPullT_Smt_Off->GetBinContent(ix + 1); //Current value
+ Double_t dCor = htmp->GetBinContent(ix + 1);
+ Double_t dRMS = htmp->GetBinError(ix + 1);
+ TH1D* hpy = fhPullT_Smt->ProjectionY(Form("%s_py%d", fhPullT_Smt->GetName(), ix), ix + 1, ix + 1);
+ if (hpy->GetEntries() > 50.) {
Int_t iBmax = hpy->GetMaximumBin();
TAxis* xaxis = hpy->GetXaxis();
Double_t dMean = xaxis->GetBinCenter(iBmax); //X-value of bin with maximal content
- Double_t dRMS = TMath::Abs(hpy->GetRMS());
+ dRMS = TMath::Abs(hpy->GetRMS());
Double_t dLim = 1.5 * dRMS;
- TFitResultPtr fRes = hpy->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
- Double_t dFMean = fRes->Parameter(1);
-
- Double_t dVal = fhPullT_Smt_Off->GetBinContent(ix + 1);
- dVal -= dFMean;
- TF1* fg = hpy->GetFunction("gaus");
- Double_t dFMeanError = fg->GetParError(1);
- LOG(info) << "Update hPullT_Smt_Off3 Ind " << ix << Form(", 0x%08x: ", fRpcAddr[ix])
- << fhPullT_Smt_Off->GetBinContent(ix + 1) << " + " << dFMean << ", Err " << dFMeanError << " -> "
- << dVal << ", Width " << dRMS << ", Chi2 " << fg->GetChisquare();
- if (dFMeanError < 0.05) { // FIXME: hardwired constant
- if (dRMS < RMSmin) dRMS = RMSmin;
- if (dRMS > fSIGT * 3.0) dRMS = fSIGT * 3.;
- if (fRpcAddr[ix] != fiBeamCounter) // don't correct beam counter time
- fhPullT_Smt_Off->SetBinContent(ix + 1, dVal);
- else
- LOG(info) << "No Off3 correction for beam counter at index " << ix;
- fhPullT_Smt_Width->SetBinContent(ix + 1, dRMS);
+ Double_t dNorm = hpy->GetBinContent(iBmax);
+ LOG(info) << "Fit3 " << hpy->GetName()
+ << Form(", %f with %f, %f, %f ", hpy->GetEntries(), dNorm, dMean, dLim);
+ if (dNorm > 10) {
+ TFitResultPtr fRes = hpy->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
+ //TF1* mgaus=new TF1("mgaus","gaus", dMean - dLim, dMean + dLim);
+ //mgaus->SetParameters(dNorm,dMean,dLim*0.5);
+ //TFitResultPtr fRes = hpy->Fit("mgaus", "SQM", "", dMean - dLim, dMean + dLim);
+ // see https://root-forum.cern.ch/t/tfitresultptr-not-valid-check/35944/4
+ if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
+ TF1* fg = hpy->GetFunction("gaus");
+ if (fg == NULL) {
+ LOG(fatal) << "No associated gaus function for " << hpy->GetName();
+ continue;
+ }
+ //Double_t dFMean = fRes->Parameter(1);
+ Double_t dFMean = fg->GetParameter(1);
+ dCor = dFMean; // update offset
+ Double_t dFMeanError = fg->GetParError(1);
+ LOG(info) << "Update hPullT_Smt_Off3 Ind " << ix << Form(", 0x%08x: ", fRpcAddr[ix])
+ << fhPullT_Smt_Off->GetBinContent(ix + 1) << " + " << dFMean << ", Err " << dFMeanError
+ << " -> " << dVal - dCor << ", Width " << dRMS << ", Chi2 " << fg->GetChisquare();
+ if (dFMeanError < 0.05) { // FIXME: hardwired constant
+ if (dRMS < RMSmin) dRMS = RMSmin;
+ if (dRMS > fSIGT * 3.0) dRMS = fSIGT * 3.;
+ }
+ }
+ else {
+ LOG(info) << " Fit of " << hpy->GetName() << " failed with " << gMinuit->fCstatu;
+ }
+ }
+ else {
+ LOG(info) << "Fit3: Too few entries for fit ofhisto " << hpy->GetName() << ": " << dNorm;
}
}
else {
- LOG(debug1) << "Update hPullT_Smt_Off " << ix << ": insufficient counts: " << hpy->GetEntries();
+ LOG(info) << "Update hPullT_Smt_Off " << ix << ": insufficient counts: " << hpy->GetEntries();
}
- }
+
+ if (fRpcAddr[ix] != fiBeamCounter) // don't correct beam counter time
+ fhPullT_Smt_Off->SetBinContent(ix + 1, dVal - dCor);
+ else
+ LOG(info) << "No Off3 correction for beam counter at index " << ix;
+
+ fhPullT_Smt_Width->SetBinContent(ix + 1, dRMS);
+ } //ix loop end
}
else {
LOG(warning) << "CbmTofFindTracks::WriteHistos: fhPullT_Smt_Off not found ";
@@ -729,6 +813,8 @@ Bool_t CbmTofFindTracks::WriteHistos()
if (fhPullX_Smt_Off != NULL) {
Double_t nx = htmp1D->GetNbinsX();
for (Int_t ix = 0; ix < nx; ix++) {
+ int iSmType = CbmTofAddress::GetSmType(fRpcAddr[ix] & DetMask);
+ if (iSmType == 8) continue; // skip pad counters
TH1D* hpy = fhPullX_Smt->ProjectionY("_py", ix + 1, ix + 1);
Double_t dVal = fhPullX_Smt_Off->GetBinContent(ix + 1);
//dVal -= htmp1D->GetBinContent(ix + 1);
@@ -774,6 +860,8 @@ Bool_t CbmTofFindTracks::WriteHistos()
if (fhPullY_Smt_Off != NULL) {
Double_t nx = htmp1D->GetNbinsX();
for (Int_t ix = 0; ix < nx; ix++) {
+ int iSmType = CbmTofAddress::GetSmType(fRpcAddr[ix] & DetMask);
+ if (iSmType == 8) continue; // skip pad counters
Double_t dVal = fhPullY_Smt_Off->GetBinContent(ix + 1);
//dVal -= htmp1D->GetBinContent(ix + 1);
// Fit gaussian
@@ -918,7 +1006,67 @@ Bool_t CbmTofFindTracks::WriteHistos()
}
} break;
- default:;
+ case 82:
+ case 81:
+ case 80: {
+ Int_t iSel = fiCorMode % 10;
+ LOG(info) << "Update time offsets with Detector Doublets " << iSel;
+ int iO[3] = {0, 5, 10}; // 0 - 1 - 2 - layers
+ switch (iSel) {
+ case 1: iO[2] = 25; break; // 0 - 1 - big modules
+ case 2:
+ iO[0] = 5;
+ iO[1] = 25;
+ iO[2] = 12;
+ break; // 1 - big - 2
+ default:;
+ }
+ const size_t N = 3;
+ double dTshift[N] = {3 * 0.};
+ for (int iLoc = 0; iLoc < 5; iLoc++) { // loop over rpcs in module
+ if (EvalDoublets(iLoc, iLoc + iO[1], iLoc + iO[2], dTshift)) { // returns vector of shifts
+ double dTMeanShift = (dTshift[0] + dTshift[1] + dTshift[2]) / 3.;
+ for (int i = 0; i < 3; i++) { // apply time shifts
+ int iStation = GetStationOfAddr(fDigiBdfPar->GetDetUId(iLoc + iO[i]));
+ if (fiStationStatus[iStation] > 0) continue; // do not modify
+ int ix = fMapRpcIdParInd[fDigiBdfPar->GetDetUId(iLoc + iO[i])]; // convert BDF to Geo counting
+ LOG(info) << "UpdateDT0 bdf ch " << iLoc + iO[i] << ", geo ch " << ix << " by " << Form("%f", dTshift[i]);
+ fhPullT_Smt_Off->SetBinContent(ix + 1, fhPullT_Smt_Off->GetBinContent(ix + 1) + dTshift[i] - dTMeanShift);
+ }
+ }
+ else {
+ iO[2]++; // try neighbor
+ if (EvalDoublets(iLoc, iLoc + iO[1], iLoc + iO[2], dTshift)) { // returns vector of shifts
+ double dTMeanShift = (dTshift[0] + dTshift[1] + dTshift[2]) / 3.;
+ for (int i = 0; i < 3; i++) {
+ int iStation = GetStationOfAddr(fDigiBdfPar->GetDetUId(iLoc + iO[i]));
+ if (fiStationStatus[iStation] > 0) continue; // do not modify
+ int ix = fMapRpcIdParInd[fDigiBdfPar->GetDetUId(iLoc + iO[i])]; // convert BDF to Geo counting
+ LOG(info) << "UpdateDT+ bdf ch " << iLoc + iO[i] << ", geo ch " << ix << " by " << dTshift[i];
+ fhPullT_Smt_Off->SetBinContent(ix + 1, fhPullT_Smt_Off->GetBinContent(ix + 1) + dTshift[i] - dTMeanShift);
+ }
+ }
+ else {
+ iO[2] -= 2; // try other neighbor
+ if (EvalDoublets(iLoc, iLoc + iO[1], iLoc + iO[2], dTshift)) { // returns vector of shifts
+ double dTMeanShift = (dTshift[0] + dTshift[1] + dTshift[2]) / 3.;
+ for (int i = 0; i < 3; i++) {
+ int iStation = GetStationOfAddr(fDigiBdfPar->GetDetUId(iLoc + iO[i]));
+ if (fiStationStatus[iStation] > 0) continue; // do not modify
+ int ix = fMapRpcIdParInd[fDigiBdfPar->GetDetUId(iLoc + iO[i])]; // convert BDF to Geo counting
+ LOG(info) << "UpdateDT- bdf ch " << iLoc + iO[i] << ", geo ch " << ix << " by " << dTshift[i];
+ fhPullT_Smt_Off->SetBinContent(ix + 1,
+ fhPullT_Smt_Off->GetBinContent(ix + 1) + dTshift[i] - dTMeanShift);
+ }
+ }
+ iO[2] += 2; // restore offset
+ }
+ iO[2]--; // restore offset
+ }
+ }
+ } break;
+
+ default: LOG(info) << "Correction mode not implemented!";
}
if (NULL != fhPullT_Smt_Off) {
@@ -976,7 +1124,7 @@ Bool_t CbmTofFindTracks::WriteHistos()
// ----- Public method Exec --------------------------------------------
void CbmTofFindTracks::Exec(Option_t* opt)
{
- if (fair::Logger::Logging(fair::Severity::debug)) { fDigiBdfPar->printParams(); }
+ //if (fair::Logger::Logging(fair::Severity::debug)) { fDigiBdfPar->printParams(); }
if (!fEventsColl) {
// fTofHitArray = (TClonesArray*)fTofHitArrayIn->Clone();
fTofHitArray = (TClonesArray*) fTofHitArrayIn;
@@ -987,6 +1135,8 @@ void CbmTofFindTracks::Exec(Option_t* opt)
Int_t iNbCalHits = 0;
fTrackArrayOut->Delete(); //Clear("C");
fTofHitArrayOut->Delete(); //Clear("C");
+ LOG(info) << "Process TS " << iTS << " with " << fEventsColl->GetEntriesFast() << " events";
+ iTS++;
for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
LOG(debug) << "Process event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kTofHit) << " hits";
@@ -1021,8 +1171,8 @@ void CbmTofFindTracks::Exec(Option_t* opt)
CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHitArray->At(iHit));
new ((*fTofHitArrayOut)[iNbCalHits++]) CbmTofHit(*tHit);
}
-
fTrackArray->Delete();
+ //fTrackArray->Clear();
}
}
}
@@ -1033,6 +1183,7 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
ResetStationsFired();
if (NULL != fTofUHitArray) fTofUHitArray->Clear("C");
if (NULL != fTrackArray) fTrackArray->Delete(); // reset
+ //if (NULL != fTrackArray) fTrackArray->Clear(); // reset
// recalibrate hits and count trackable hits
for (Int_t iHit = 0; iHit < fTofHitArray->GetEntriesFast(); iHit++) {
@@ -1089,45 +1240,37 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
}
// tune positions and times
- Double_t dTcor = 0.;
if ((iDetId & DetMask) != fiBeamCounter) { // do not modify diamond position
Int_t iRpcInd = fMapRpcIdParInd[iDetId];
pHit->SetTime(pHit->GetTime() + fvToff[iRpcInd]);
pHit->SetX(pHit->GetX() + fvXoff[iRpcInd]);
pHit->SetY(pHit->GetY() + fvYoff[iRpcInd]);
pHit->SetZ(pHit->GetZ() + fvZoff[iRpcInd]);
- /*
- if (fhPullT_Smt_Off != NULL) {
- dTcor = (Double_t) fhPullT_Smt_Off->GetBinContent(iRpcInd + 1);
- pHit->SetTime(pHit->GetTime() + dTcor);
- }
- if (fhPullX_Smt_Off != NULL) {
- Double_t dXcor = (Double_t) fhPullX_Smt_Off->GetBinContent(iRpcInd + 1);
- pHit->SetX(pHit->GetX() + dXcor);
- }
- if (fhPullY_Smt_Off != NULL) {
- Double_t dYcor = (Double_t) fhPullY_Smt_Off->GetBinContent(iRpcInd + 1);
- pHit->SetY(pHit->GetY() + dYcor);
- }
- if (fhPullZ_Smt_Off != NULL) {
- Double_t dZcor = (Double_t) fhPullZ_Smt_Off->GetBinContent(iRpcInd + 1);
- pHit->SetZ(pHit->GetZ() + dZcor);
- }
- */
}
Int_t iSt = GetStationOfAddr(iDetId);
- MarkStationFired(iSt);
+ if (iSt >= (Int_t) fStationHMul.size()) {
+ LOG(fatal) << Form("Invalid station # %d for detId 0x%08x, TSR %d%d%d", iSt, iDetId,
+ CbmTofAddress::GetSmType(iDetId), CbmTofAddress::GetSmId(iDetId),
+ CbmTofAddress::GetRpcId(iDetId));
+ }
+
+ if ((Int_t) fMapStationRpcId.size() > fNTofStations) {
+ PrintSetup();
+ LOG(fatal) << "Invalid NTofStations " << fNTofStations << ", " << fMapStationRpcId.size();
+ }
- LOG(debug) << Form("Exec found Hit %02d, addr 0x%08x, sta %02d, %02d, HM "
+ LOG(debug) << Form("Exec found 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)",
- iHit, pHit->GetAddress(), GetStationOfAddr(iDetId), fDigiBdfPar->GetTrackingStation(pHit),
+ "%6.2f(%3.2f) Z %6.2f(%3.2f) T %6.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(), dTcor);
+ pHit->GetZ(), pHit->GetDz(), pHit->GetTime(), pHit->GetTimeError());
+ MarkStationFired(iSt);
}
- LOG(debug) << Form("CbmTofFindTracks::Exec NStationsFired %d > %d Min ?, NbStations %d", GetNStationsFired(),
+ LOG(debug) << Form("CbmTofFindTracks::Exec NStationsFired %d >= %d Min ?, NbStations %d", GetNStationsFired(),
GetMinNofHits(), fDigiBdfPar->GetNbTrackingStations());
if (GetNStationsFired() < GetMinNofHits()) {
@@ -1153,7 +1296,7 @@ void CbmTofFindTracks::ExecFind(Option_t* /*opt*/, CbmEvent* tEvent)
FindVertex();
- FillHistograms(tEvent);
+ if (fbDoHistos) FillHistograms(tEvent);
}
FillUHits(); // put unused hits into TClonesArray
@@ -1166,7 +1309,7 @@ void CbmTofFindTracks::Finish()
{
if (fiEvent < 1000) return; // preserve calibration histos in event display
if (fiCalOpt > 0) fTofCalibrator->UpdateCalHist(fiCalOpt);
- WriteHistos();
+ if (fbDoHistos) WriteHistos();
LOG(info) << Form(" CbmTofFindTracks::Finished ");
}
@@ -1196,10 +1339,10 @@ void CbmTofFindTracks::CreateHistograms()
((CbmTofTrackFinderNN*) fFinder)->GetChiMaxAccept());
fhTrackingTimeNhits =
- new TH2F(Form("hTrackingTimeNhits"), Form("Tracking Time; NHits; #Deltat (s)"), 100, 0, 200, 50, 0, 0.1);
+ new TH2F(Form("hTrackingTimeNhits"), Form("Tracking Time; NHits; #Deltat (s)"), 200, 0, 200, 50, 0, 0.2);
fhTrklMulNhits =
- new TH2F(Form("hTrklMulNhits"), Form("Tracklet Multiplicity; NHits; NTracklet"), 150, 0, 150, 30, 0, 30);
+ new TH2F(Form("hTrklMulNhits"), Form("Tracklet Multiplicity; NHits; NTracklet"), 200, 0, 200, 30, 0, 30);
fhTrklMulMaxMM =
new TH2F(Form("hTrklMulMaxMax-1"), Form("Tracklet Multiplicity; TMulMax; TMulMax-1"), 10, 0, 10, 10, 0, 10);
@@ -1213,11 +1356,22 @@ void CbmTofFindTracks::CreateHistograms()
fhTrklZ0yHMul = new TH2F(Form("hTrklZ0yHMul"), Form("Tracklet Z0y vs. Hit - Multiplicity; HMul_{Tracklet}; Z0y"), 8,
2, 10, 100, -300, 300);
+ double dTxLim = ((CbmTofTrackFinderNN*) fFinder)->GetTxLIM();
+ double dTyLim = ((CbmTofTrackFinderNN*) fFinder)->GetTyLIM();
+ double dTxMean = ((CbmTofTrackFinderNN*) fFinder)->GetTxMean();
+ double dTyMean = ((CbmTofTrackFinderNN*) fFinder)->GetTyMean();
+ double dTxMin = dTxMean - dTxLim;
+ double dTxMax = dTxMean + dTxLim;
+ double dTyMin = dTyMean - dTyLim;
+ double dTyMax = dTyMean + dTyLim;
+
fhTrklTxHMul = new TH2F(Form("hTrklTxHMul"), Form("Tracklet Tx vs. Hit - Multiplicity; HMul_{Tracklet}; Tx"), 8, 2,
- 10, 100, -0.65, 0.65);
+ 10, 100, dTxMin, dTxMax);
+ fhTrklTyHMul = new TH2F(Form("hTrklTyHMul"), Form("Tracklet Ty vs. Hit - Multiplicity; HMul_{Tracklet}; Ty"), 8, 2,
+ 10, 100, dTyMin, dTyMax);
+
+ fhTrklTyTx = new TH2F(Form("hTrklTyTx"), Form("Tracklet Ty vs. Tx; Tx; Ty"), 50, dTxMin, dTxMax, 50, dTyMin, dTyMax);
- fhTrklTyHMul = new TH2F(Form("hTrklTyHMul"), Form("Tracklet Ty vs. Hit - Multiplicity; HMul_{Tracklet}; Ty"), 8, 2,
- 10, 100, -0.65, 0.65);
Double_t TTMAX = 0.2;
fhTrklTtHMul = new TH2F(Form("hTrklTtHMul"), Form("Tracklet Tt vs. Hit - Multiplicity; HMul_{Tracklet}; Tt"), 8, 2,
10, 100, -TTMAX, TTMAX);
@@ -1263,7 +1417,7 @@ void CbmTofFindTracks::CreateHistograms()
fhTOff_HMul2 = new TH2F(Form("hTOff_HMul2"), Form("Tracklet TOff(HMul2); RpcInd ; TOff (ns)"), nSmt, 0, nSmt, 500,
-fT0MAX, fT0MAX);
- Double_t DTTMAX = 0.09;
+ Double_t DTTMAX = 0.11;
fhDeltaTt_Smt = new TH2F(Form("hDeltaTt_Smt"), Form("Tracklet DeltaTt; RpcInd ; #DeltaTt (ns/cm)"), nSmt, 0, nSmt,
100, -DTTMAX, DTTMAX);
@@ -1272,10 +1426,12 @@ void CbmTofFindTracks::CreateHistograms()
vhPullZ.resize(fNTofStations);
vhPullT.resize(fNTofStations);
vhPullTB.resize(fNTofStations);
+ vhTrefRms.resize(fNTofStations);
vhResidualTBWalk.resize(fNTofStations);
vhResidualYWalk.resize(fNTofStations);
vhXY_AllTracks.resize(fNTofStations);
vhXY_AllStations.resize(fNTofStations);
+ vhXY_AllFitStations.resize(fNTofStations);
vhXY_MissedStation.resize(fNTofStations);
vhXY_DX.resize(fNTofStations);
vhXY_DY.resize(fNTofStations);
@@ -1296,6 +1452,7 @@ void CbmTofFindTracks::CreateHistograms()
new TH1F(Form("hPullT_Station_%d", iSt), Form("hResiT_Station_%d; #DeltaT (ns)", iSt), 59, -fT0MAX, fT0MAX);
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.);
const Double_t TOTmax = 50.;
vhResidualTBWalk[iSt] =
new TH2F(Form("hResidualTBWalk_Station_%d", iSt), Form("hResidualTBWalk_Station_%d; #DeltaT (ns)", iSt), TOTmax,
@@ -1307,7 +1464,7 @@ void CbmTofFindTracks::CreateHistograms()
Int_t Nbins = 32.;
CbmTofCell* fChannelInfo = fDigiPar->GetCell(fMapStationRpcId[iSt]);
if (NULL == fChannelInfo) {
- LOG(fatal) << "Geometry for station " << iSt << ", Rpc " << fMapStationRpcId[iSt] << " not defined ";
+ LOG(fatal) << "Geometry for station " << iSt << Form(", Rpc %0x08x ", fMapStationRpcId[iSt]) << " not defined ";
return;
}
Int_t NbinsX = 2 * XSIZ / fChannelInfo->GetSizex();
@@ -1315,6 +1472,9 @@ void CbmTofFindTracks::CreateHistograms()
-XSIZ, XSIZ, Nbins, -XSIZ, XSIZ);
vhXY_AllStations[iSt] = new TH2F(Form("hXY_AllStations_%d", iSt), Form("hXY_AllStations_%d; x(cm); y (cm)", iSt),
NbinsX, -XSIZ, XSIZ, Nbins, -XSIZ, XSIZ);
+ vhXY_AllFitStations[iSt] =
+ new TH2F(Form("hXY_AllFitStations_%d", iSt), Form("hXY_AllFitStations_%d; x(cm); y (cm)", iSt), NbinsX, -XSIZ,
+ XSIZ, Nbins, -XSIZ, XSIZ);
vhXY_MissedStation[iSt] =
new TH2F(Form("hXY_MissedStation_%d", iSt), Form("hXY_MissedStation_%d; x(cm); y (cm)", iSt), NbinsX, -XSIZ,
XSIZ, Nbins, -XSIZ, XSIZ);
@@ -1341,7 +1501,7 @@ void CbmTofFindTracks::CreateHistograms()
}
- // vertex histrograms
+ // vertex histograms
Double_t NNORM = 40.;
fhVTXNorm = new TH1F(Form("hVTXNorm"), Form("Vertex Normalisation; #_{TrackletHits}"), NNORM, 0, NNORM);
fhVTX_XY0 =
@@ -1408,6 +1568,8 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
}
}
if (fiBeamCounter != -1 && NULL == pRefHit) return;
+ //LOG(debug)<<"FillH_start";
+ //PrintSetup();
std::vector<Int_t> HMul;
HMul.resize(fNTofStations + 1);
@@ -1419,13 +1581,15 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
for (Int_t iTrk = 0; iTrk < fTrackArray->GetEntriesFast(); iTrk++) {
CbmTofTracklet* pTrk = (CbmTofTracklet*) fTrackArray->At(iTrk);
if (NULL == pTrk) continue;
- if (pTrk->GetNofHits() > fNTofStations) {
+ if (pTrk->GetNofHits() > fNTofStations + 1) { // +1 for virtual T0
LOG(error) << "CbmTofFindTracks::FillHistograms: more hits (" << pTrk->GetNofHits() << ") than stations ("
<< fNTofStations << ")";
continue;
}
- //if (pTrk->GetTrackTy()>0.) continue; // TEMPORARY!!!!
+ //if (pTrk->GetTrackTy()>0.) continue; // Debugging!!!!
+ //LOG(debug)<<"FillH_HMul "<<iTrk<<", "<<pTrk->GetNofHits();
+ //PrintSetup();
HMul[pTrk->GetNofHits()]++;
@@ -1441,6 +1605,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
Int_t iDetId = pTrk->GetTofDetIndex(iH); // DetId of iH. Hit
CbmTofHit* pHit = pTrk->GetTofHitPointer(iH);
Int_t iSt = GetStationOfAddr(iDetId); // Station of 1. Hit
+ if (iSt < 0) continue;
Double_t dTOff = pHit->GetTime() - pHit->GetR() * fTtTarg - dTRef0;
LOG(debug1) << Form("<D> CbmTofFindTracks::FillHistograms: iDetId1 "
"0x%08x, iST1 = %d with dTOff %f at RpcInd %d",
@@ -1449,15 +1614,16 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
} // loop over tracklets' hits
}
- if (pTrk->GetNofHits() > fMinNofHits) { // for further analysis request at least 3 matched hits
+ if (pTrk->GetNofHits() >= fMinNofHits) { // for further analysis request at least 3 matched hits
iTMul++;
fhTrklChi2->Fill(pTrk->GetNofHits(), pTrk->GetChiSq());
- if (fiCalOpt > 0) fTofCalibrator->FillCalHist(pTrk, fiCalOpt, tEvent);
+ if (fiCalOpt > 0)
+ if (pTrk->GetTt() > fdTtMin) fTofCalibrator->FillCalHist(pTrk, fiCalOpt, tEvent);
CbmTofTrackletParam* tPar = pTrk->GetTrackParameter();
Double_t dTt = pTrk->GetTt();
- LOG(debug) << Form("Trk %d info: Lz=%6.2f Z0x=%6.2f Z0y=%6.2f Tt=%6.4f", iTrk, tPar->GetLz(), pTrk->GetZ0x(),
+ LOG(debug) << Form("Trk %d info: Lz=%6.2f Z0x=%6.2f Z0y=%6.2f Tt=%6.4f, ", iTrk, tPar->GetLz(), pTrk->GetZ0x(),
pTrk->GetZ0y(), dTt)
<< tPar->ToString();
@@ -1465,8 +1631,20 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
fhTrklZ0yHMul->Fill(pTrk->GetNofHits(), pTrk->GetFitY(0.));
fhTrklTxHMul->Fill(pTrk->GetNofHits(), tPar->GetTx());
fhTrklTyHMul->Fill(pTrk->GetNofHits(), tPar->GetTy());
+ fhTrklTyTx->Fill(tPar->GetTx(), tPar->GetTy());
fhTrklTtHMul->Fill(pTrk->GetNofHits(), dTt);
+ if (kFALSE) { //kFALSE && fdTtMin < -10.) { // print event number for inspection in event display
+ //Double_t zPosTar = 160.;
+ //if (TMath::Abs(tPar->GetTx())<0.2 && TMath::Abs(tPar->GetTy())<0.1) {
+ //if(TMath::Abs(pTrk->GetFitX(0.))<5 && TMath::Abs(pTrk->GetFitY(0.))<5) {
+ //if ((TMath::Abs(pTrk->GetFitX(zPosTar) - 20.) < 10 && TMath::Abs(pTrk->GetFitY(zPosTar) - 15.) < 5)
+ if (pTrk->GetNofHits() > 8) {
+ LOG(info) << Form("Found tracklet candidate (x0,y0)=(%5.2f,%5.2f), Hmul %d, v %6.2f in event %d ",
+ pTrk->GetFitX(0.), pTrk->GetFitY(0.), pTrk->GetNofHits(), 1. / dTt, fiEvent)
+ << " within " << fTofHitArray->GetEntriesFast() << " hits ";
+ }
+ }
switch (GetNReqStations() - pTrk->GetNofHits()) {
case 0: // max hit number
fhTrklXY0_0->Fill(pTrk->GetFitX(0.), pTrk->GetFitY(0.));
@@ -1481,7 +1659,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
if (pTrk->GetR0() > fdR0Lim) continue;
}
- if (dTt > 0.)
+ if (dTt > fdTtMin)
for (Int_t iSt = 0; iSt < fNTofStations; iSt++) {
Int_t iH = pTrk->GetStationHitIndex(fMapStationRpcId[iSt]); // Station Hit index
if (iH < 0) continue; // Station not part of tracklet
@@ -1498,8 +1676,10 @@ 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 dDTB = fTrackletTools->GetTdif(pTrk, fMapStationRpcId[iSt],
- pHit); // ignore pHit in calc of reference
+ Double_t dTexp = fTrackletTools->GetTexpected(pTrk, fMapStationRpcId[iSt], pHit);
+ Double_t dDTB = pHit->GetTime() - dTexp;
+ // fTrackletTools->GetTdif(pTrk, fMapStationRpcId[iSt],pHit); // ignore pHit in calc of reference
+ Double_t dTexpErr = fTrackletTools->GetTexpectedError(pTrk, fMapStationRpcId[iSt], pHit, dTexp);
Double_t dTOT = pHit->GetCh() / 10.; // misuse of channel field
Double_t dZZ = pHit->GetZ() - tPar->GetZy(pHit->GetY());
@@ -1513,6 +1693,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
vhPullZ[iSt]->Fill(dZZ);
vhPullT[iSt]->Fill(dDT);
vhPullTB[iSt]->Fill(dDTB);
+ vhTrefRms[iSt]->Fill(dTexpErr);
vhResidualTBWalk[iSt]->Fill(dTOT, dDTB);
vhResidualYWalk[iSt]->Fill(dTOT, dDY);
@@ -1530,20 +1711,26 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
fhDeltaTt_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dDeltaTt);
//XXX use BRef as Referenz!!!
if (pRefHit != NULL) {
- Double_t dTOff = dDeltaTt * //pHit->GetR();
- TMath::Sqrt(TMath::Power(pHit->GetX() - pRefHit->GetX(), 2)
- + TMath::Power(pHit->GetY() - pRefHit->GetY(), 2)
- + TMath::Power(pHit->GetZ() - pRefHit->GetZ(), 2))
- * TMath::Sign(1, pHit->GetZ() - pRefHit->GetZ());
- fhTOff_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dTOff);
+ if (pHit != pRefHit) {
+ Double_t dTOff = dDeltaTt * //pHit->GetR();
+ TMath::Sqrt(TMath::Power(pHit->GetX() - pRefHit->GetX(), 2)
+ + TMath::Power(pHit->GetY() - pRefHit->GetY(), 2)
+ + TMath::Power(pHit->GetZ() - pRefHit->GetZ(), 2))
+ * TMath::Sign(1, pHit->GetZ() - pRefHit->GetZ());
+ fhTOff_Smt->Fill((Double_t) fMapRpcIdParInd[fMapStationRpcId[iSt]], dTOff);
+ }
}
- }
+ } // station loop end
// extrapolation of tracklet to vertex @ z=0
// FairTrackParam paramExtr;
// fFitter->Extrapolate(pTrk->GetParamFirst(),0.,¶mExtr);
- } // condition on NofHits>2 end
-
+ } // condition on NofHits>2 end
+ Double_t dTt = fTrackletTools->FitTt(pTrk, -1); // restore full fit
+ if (dTt == 0) {
+ LOG(error) << "Invalid inverse velocity ";
+ continue;
+ }
// monitoring of tracklet hits with selected velocities from reference counters
if (TMath::Abs(pTrk->GetRefVel((UInt_t)(fNReqStations - 1)) - fdRefVelMean) < fdRefDVel) {
fhTrklVelHMul->Fill(pTrk->GetNofHits(), 1. / pTrk->GetTt());
@@ -1571,7 +1758,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
for (Int_t iSt = 0; iSt < fNTofStations; iSt++) {
Int_t iH = pTrk->GetStationHitIndex(fMapStationRpcId[iSt]); // Station Hit index
if (iH < 0) {
- LOG(debug) << " Incomplete Tracklet, skip station " << iSt;
+ LOG(debug1) << " Incomplete Tracklet, skip station " << iSt;
continue; // Station not part of tracklet
}
CbmTofHit* pHit = (CbmTofHit*) fTofHitArray->At(iH);
@@ -1606,8 +1793,17 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
vhXY_TOT[iSt]->Fill(hitpos_local[0], hitpos_local[1], dTOT);
vhXY_CSZ[iSt]->Fill(hitpos_local[0], hitpos_local[1], dCSZ);
+ // Refit tracklet without Dut info
+ fTrackletTools->Line3DFit(pTrk, iChId & DetMask);
+ hitpos[0] = pTrk->GetFitX(pHit->GetZ());
+ hitpos[1] = pTrk->GetFitY(pHit->GetZ());
+ hitpos[2] = pHit->GetZ();
+ gGeoManager->MasterToLocal(hitpos, hitpos_local);
+ vhXY_AllFitStations[iSt]->Fill(hitpos_local[0], hitpos_local[1]);
+ fTrackletTools->Line3DFit(pTrk, -1);
+
// debugging consistency of geometry transformations ....
- if (fair::Logger::Logging(fair::Severity::debug)) {
+ if (fair::Logger::Logging(fair::Severity::debug3)) {
if (iSt == fNReqStations - 1) { // treat as if not found
Int_t iAddr = fMapStationRpcId[iSt];
CbmTofCell* fChannelInfoD = fDigiPar->GetCell(iAddr);
@@ -1719,11 +1915,12 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
fChannelInfo->GetSizex());
}
zPosMiss = fChannelInfoMiss->GetZ() + (Double_t) fhPullZ_Smt_Off->GetBinContent(iRpcInd + 1);
-
+ /*
LOG(debug) << Form("Geo consistency check 0x%08x at St%d, z=%7.2f,%7.2f: "
"iChTrafo %d, Miss %d , xloc %6.2f, dx %4.2f",
iAddr, iSt, zPos, zPosMiss, iChTrafo, iChMiss, hitpos_local[0],
fChannelInfo->GetSizex());
+ */
fChannelInfo = fChannelInfoMiss;
iAddr = iAddrMiss;
}
@@ -1764,7 +1961,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
} // velocity selection end
} // loop over tracklets end
- if (HMul.size() > 3) fhTrklMulMaxMM->Fill(HMul[fNTofStations], HMul[fNTofStations - 1]);
+ if (HMul.size() > 3) { fhTrklMulMaxMM->Fill(HMul[fNReqStations], HMul[fNReqStations - 1]); }
if (HMul.size() > 5) fhTrklMul3D->Fill(HMul[fNTofStations], HMul[fNTofStations - 1], HMul[fNTofStations - 2]);
fhTrklMulNhits->Fill(fTofHitArray->GetEntriesFast(), iTMul);
fhTrackingTimeNhits->Fill(fTofHitArray->GetEntriesFast(), fdTrackingTime);
@@ -1783,7 +1980,7 @@ void CbmTofFindTracks::FillHistograms(CbmEvent* tEvent)
pTrk->PrintInfo();
}
}
- if (1)
+ if (kFALSE) // print event info for special events
if (fTrackArray->GetEntriesFast() > 25) { // temporary
LOG(info) << "Found high track multiplicity of " << fTrackArray->GetEntriesFast() << " in event " << fiEvent
<< " from " << fTofHitArray->GetEntriesFast() << " hits ";
@@ -1869,6 +2066,8 @@ void CbmTofFindTracks::SetStation(Int_t iVal, Int_t iModType, Int_t iModId, Int_
if (NULL != fDigiBdfPar) iCenterCh = TMath::Floor((fDigiBdfPar->GetNbChan(iModType, iRpcId) - 1) / 2);
Int_t iAddr = CbmTofAddress::GetUniqueAddress(iModId, iRpcId, iCenterCh, 0, iModType);
fMapStationRpcId[iVal] = iAddr;
+ LOG(info) << "SetStation: added " << iVal
+ << Form(" TSRC %d%d%d%02d, addr 0x%08x ", iModType, iModId, iRpcId, iCenterCh, iAddr);
}
void CbmTofFindTracks::SetBeamCounter(Int_t iModType, Int_t iModId, Int_t iRpcId)
@@ -1879,9 +2078,21 @@ void CbmTofFindTracks::SetBeamCounter(Int_t iModType, Int_t iModId, Int_t iRpcId
Int_t CbmTofFindTracks::GetStationOfAddr(Int_t iAddr)
{
std::map<Int_t, Int_t>::iterator it;
- for (it = fMapStationRpcId.begin(); it != fMapStationRpcId.end(); it++)
+ Int_t iNSt = 0;
+ Int_t iSt = -1;
+ for (it = fMapStationRpcId.begin(); it != fMapStationRpcId.end(); it++) {
//std::map <Int_t, Int_t>::iterator it = fMapStationRpcId.find(iAddr);
+ /*
+ Int_t iDetId=it->second;
+ LOG(debug)<<Form("GetStationOfAddr 0x%08x, TSR %d%d%d: St %d, TSR %d%d%d ", iAddr,
+ CbmTofAddress::GetSmType(iAddr),CbmTofAddress::GetSmId(iAddr),CbmTofAddress::GetRpcId(iAddr),
+ it->first, CbmTofAddress::GetSmType(iDetId),CbmTofAddress::GetSmId(iDetId),CbmTofAddress::GetRpcId(iDetId));
+ */
+ if (iNSt > fNTofStations) LOG(fatal) << "Invalid NSt " << iNSt;
+ iNSt++;
if (it->second == iAddr) break;
+ }
+ if (it != fMapStationRpcId.end()) iSt = it->first;
/*
if(it->first == fMapStationRpcId.size())
{
@@ -1890,14 +2101,19 @@ Int_t CbmTofFindTracks::GetStationOfAddr(Int_t iAddr)
;
}
*/
- return it->first;
+ return iSt;
}
void CbmTofFindTracks::PrintSetup()
{
for (std::map<Int_t, Int_t>::iterator it = fMapStationRpcId.begin(); it != fMapStationRpcId.end(); it++) {
- cout << " <I> Tracking station " << it->first << " contains RpcId " << Form("0x%08x", it->second) << endl;
+ Int_t iDetId = it->second;
+ LOG(debug) << " <I> Tracking station " << it->first << " contains RpcId "
+ << Form("0x%08x, TSR %d%d%d", it->second, CbmTofAddress::GetSmType(iDetId),
+ CbmTofAddress::GetSmId(iDetId), CbmTofAddress::GetRpcId(iDetId));
}
+ if ((Int_t) fMapStationRpcId.size() > fNTofStations)
+ LOG(fatal) << "Invalid NTofStations " << fNTofStations << ", " << fMapStationRpcId.size();
}
Double_t CbmTofFindTracks::GetTOff(Int_t iAddr)
@@ -1907,32 +2123,28 @@ Double_t CbmTofFindTracks::GetTOff(Int_t iAddr)
return (Double_t) fhPullT_Smt_Off->GetBinContent(fMapRpcIdParInd[iAddr] + 1);
}
-Double_t CbmTofFindTracks::GetSigT(Int_t iAddr)
-{
- return (Double_t) fhPullT_Smt_Width->GetBinContent(fMapRpcIdParInd[iAddr] + 1);
-}
+Double_t CbmTofFindTracks::GetSigT(Int_t iAddr) { return fvTsig[GetStationOfAddr(iAddr)]; }
-Double_t CbmTofFindTracks::GetSigX(Int_t iAddr)
-{
- return (Double_t) fhPullX_Smt_Width->GetBinContent(fMapRpcIdParInd[iAddr] + 1);
-}
+Double_t CbmTofFindTracks::GetSigX(Int_t iAddr) { return fvXsig[GetStationOfAddr(iAddr)]; }
-Double_t CbmTofFindTracks::GetSigY(Int_t iAddr)
-{
- return (Double_t) fhPullY_Smt_Width->GetBinContent(fMapRpcIdParInd[iAddr] + 1);
-}
+Double_t CbmTofFindTracks::GetSigY(Int_t iAddr) { return fvYsig[GetStationOfAddr(iAddr)]; }
-Double_t CbmTofFindTracks::GetSigZ(Int_t iAddr)
-{
- return (Double_t) fhPullZ_Smt_Width->GetBinContent(fMapRpcIdParInd[iAddr] + 1);
-}
+Double_t CbmTofFindTracks::GetSigZ(Int_t iAddr) { return fvZsig[GetStationOfAddr(iAddr)]; }
+
+Double_t CbmTofFindTracks::GetStationSigT(Int_t iSt) { return fvTsig[iSt]; }
+
+Double_t CbmTofFindTracks::GetStationSigX(Int_t iSt) { return fvXsig[iSt]; }
+
+Double_t CbmTofFindTracks::GetStationSigY(Int_t iSt) { return fvYsig[iSt]; }
+
+Double_t CbmTofFindTracks::GetStationSigZ(Int_t iSt) { return fvZsig[iSt]; }
Int_t CbmTofFindTracks::GetNStationsFired()
{
Int_t iNSt = 0;
for (Int_t iSt = 0; iSt < fNTofStations; iSt++) {
if (fStationHMul[iSt] > 0 && fStationHMul[iSt] < fiStationMaxHMul) iNSt++;
- LOG(debug2) << Form("Station %d, Addr 0x%08x, HMul %d, Max %d, Sum %d", iSt, GetAddrOfStation(iSt),
+ LOG(debug1) << Form("Station %d, Addr 0x%08x, HMul %d, Max %d, Sum %d", iSt, GetAddrOfStation(iSt),
fStationHMul[iSt], fiStationMaxHMul, iNSt);
}
return iNSt;
@@ -1985,8 +2197,7 @@ Bool_t CbmTofFindTracks::CheckHit2Track(CbmTofHit* pHit)
void CbmTofFindTracks::CheckMaxHMul()
{
- fInspectEvent = kTRUE;
-
+ //fInspectEvent = kTRUE;
for (Int_t iSt = 0; iSt < fNTofStations; iSt++) {
if (fStationHMul[iSt] > fiStationMaxHMul) { fInspectEvent = kFALSE; }
else {
@@ -2003,3 +2214,75 @@ void CbmTofFindTracks::PrintMapRpcIdParInd()
it++;
}
}
+
+void CbmTofFindTracks::MarkStationFired(Int_t iSt) { fStationHMul[iSt]++; }
+
+bool CbmTofFindTracks::EvalDoublets(int iI0, int iI1, int iI2, double* dTshift)
+{
+ //LOG(info)<<"Evaluate time offsets from 3 doublets in triplet "<<iI0<<" "<<iI1<<" "<<iI2;
+ //for (int i=0; i<3; i++) dTshift[i]=(double)i;
+ //LOG(info)<<"Initial Tshift " <<dTshift[0]<<" "<<dTshift[1]<<" "<<dTshift[2];
+ //return kTRUE;
+
+ const double c = 30.; //speed of light in cm/ns
+ const double dNTrks = 100.;
+ double A[3] = {3 * 0.};
+ double D[3] = {3 * 0.};
+ TString hnameDT[3] = {Form("hDoubletDt_%02d%02d", iI0, iI1), Form("hDoubletDt_%02d%02d", iI0, iI2),
+ Form("hDoubletDt_%02d%02d", iI1, iI2)};
+ TString hnameDD[3] = {Form("hDoubletDd_%02d%02d", iI0, iI1), Form("hDoubletDd_%02d%02d", iI0, iI2),
+ Form("hDoubletDd_%02d%02d", iI1, iI2)};
+ TString hnameV[3] = {Form("hDoubletV_%02d%02d", iI0, iI1), Form("hDoubletV_%02d%02d", iI0, iI2),
+ Form("hDoubletV_%02d%02d", iI1, iI2)};
+ TH1* hDT[3];
+ TH1* hDD[3];
+ TH1* hV[3];
+ for (int i = 0; i < 3; i++) {
+ hDT[i] = (TH1*) gROOT->FindObjectAny(hnameDT[i]);
+ if (NULL == hDT[i]) {
+ LOG(warn) << "CalHisto " << hnameDT[i] << " not existing";
+ return kFALSE;
+ }
+ hDD[i] = (TH1*) gROOT->FindObjectAny(hnameDD[i]);
+ if (NULL == hDD[i]) {
+ LOG(warn) << "CalHisto " << hnameDD[i] << " not existing";
+ return kFALSE;
+ }
+ hV[i] = (TH1*) gROOT->FindObjectAny(hnameV[i]);
+ if (NULL == hV[i]) {
+ LOG(warn) << "CalHisto " << hnameV[i] << " not existing";
+ return kFALSE;
+ }
+ if (hDT[i]->Integral() < dNTrks) {
+ LOG(warn) << "Too few entries for triplet " << iI0 << " " << iI1 << " " << iI2 << " in " << hnameV[i] << ": "
+ << hV[i]->Integral();
+ return kFALSE; // return without shifts
+ }
+ //D[i]=hDT[i]->GetMean()*hV[i]->GetMean();
+ D[i] = hDD[i]->GetMean();
+ A[i] = (D[i] - c * hDT[i]->GetMean()) / c;
+ }
+ dTshift[0] = 0.5 * (A[0] + A[1] - A[2]);
+ dTshift[1] = 0.5 * (-A[0] + A[1] - A[2]);
+ dTshift[2] = 0.5 * (A[0] + A[1] + A[2]);
+ LOG(info) << "Final Tshift for " << iI0 << " " << iI1 << " " << iI2 << ": " << dTshift[0] << " " << dTshift[1] << " "
+ << dTshift[2];
+ return kTRUE;
+}
+
+void CbmTofFindTracks::SetStationStatus(int iStation, int iStatus)
+{
+ if ((int) fiStationStatus.size() < iStation) {
+ fiStationStatus.resize(iStation + 1);
+ for (int i = 0; i < iStation + 1; i++)
+ fiStationStatus[i] = 0;
+ LOG(info) << "StationStatus initialized for " << fiStationStatus.size() << " stations";
+ }
+ fiStationStatus[iStation] = iStatus;
+}
+int CbmTofFindTracks::GetStationStatus(int iStation)
+{
+ if (iStation < (int) fiStationStatus.size()) return fiStationStatus[iStation];
+ else
+ return 0;
+}
diff --git a/reco/detectors/tof/CbmTofFindTracks.h b/reco/detectors/tof/CbmTofFindTracks.h
index 363b867d0cad6f1d0e7cefcddf261f5fc04fa2d6..4633ac55e4c5e01a6a0d8e185a7082e40f1b0f6a 100644
--- a/reco/detectors/tof/CbmTofFindTracks.h
+++ b/reco/detectors/tof/CbmTofFindTracks.h
@@ -55,7 +55,7 @@ public:
CbmTofFindTracks();
- /** Standard constructor
+ /** Standard constructor
**
*@param name Name of class
*@param title Task title
@@ -118,6 +118,8 @@ public:
void PrintSetup();
inline void SetR0Lim(Double_t dVal) { fdR0Lim = dVal; }
+ inline Double_t GetR0Lim() { return fdR0Lim; }
+ inline void SetTtMin(Double_t dVal) { fdTtMin = dVal; }
inline Int_t GetAddrOfStation(Int_t iVal) { return fMapStationRpcId[iVal]; }
inline Int_t GetDetIndSize() { return fMapRpcIdParInd.size(); }
Int_t GetStationOfAddr(Int_t iAddr);
@@ -129,6 +131,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 GetTOffScal() const { return fdTOffScal; }
inline Double_t GetSigT() const { return fSIGT; }
inline Double_t GetSigX() const { return fSIGX; }
@@ -142,6 +145,11 @@ public:
Double_t GetSigZ(Int_t iAddr);
Double_t GetTOff(Int_t iAddr);
+ Double_t GetStationSigT(Int_t iSt);
+ Double_t GetStationSigX(Int_t iSt);
+ Double_t GetStationSigY(Int_t iSt);
+ Double_t GetStationSigZ(Int_t iSt);
+
inline void SetSIGT(Double_t dval) { fSIGT = dval; }
inline void SetSIGX(Double_t dval) { fSIGX = dval; }
inline void SetSIGY(Double_t dval) { fSIGY = dval; }
@@ -154,17 +162,21 @@ 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 SetTOffScal(Double_t val) { fdTOffScal = val; }
inline void SetT0MAX(Double_t val) { fT0MAX = val; }
inline void SetStationMaxHMul(Int_t ival) { fiStationMaxHMul = ival; }
- inline void MarkStationFired(Int_t iSt) { fStationHMul[iSt]++; }
+ void MarkStationFired(Int_t iSt);
Int_t GetNStationsFired();
void ResetStationsFired();
+ void SetStationStatus(int iStation, int iStatus);
+ int GetStationStatus(int iStation);
inline void SetBeamMomentumLab(Double_t dval) { fdBeamMomentumLab = dval; }
inline void SetRemoveSignalPropagationTime(Bool_t bval) { fbRemoveSignalPropagationTime = bval; }
inline void SetBeamMaxHMul(Int_t ival) { fiBeamMaxHMul = ival; }
inline void SetCalOpt(Int_t ival) { fiCalOpt = ival; }
+ inline void SetNoHistos() { fbDoHistos = kFALSE; }
inline Double_t GetVertexT() const { return fVTX_T; }
inline Double_t GetVertexX() const { return fVTX_X; }
@@ -177,6 +189,7 @@ public:
else
return fTofHitIndexArray[iHit];
}
+ bool EvalDoublets(int iI0, int iI1, int iI2, double* dTshift);
private:
static CbmTofFindTracks* fInstance;
@@ -211,6 +224,11 @@ private:
std::vector<Double_t> fvYoff; // station correction parameter
std::vector<Double_t> fvZoff; // station correction parameter
+ std::vector<Double_t> fvTsig; // station resolution parameter
+ std::vector<Double_t> fvXsig; // station resolution parameter
+ std::vector<Double_t> fvYsig; // station resolution parameter
+ std::vector<Double_t> fvZsig; // station resolution parameter
+
CbmTofFindTracks(const CbmTofFindTracks&);
CbmTofFindTracks& operator=(const CbmTofFindTracks&);
@@ -232,6 +250,7 @@ private:
TH2* fhTrklZ0yHMul;
TH2* fhTrklTxHMul;
TH2* fhTrklTyHMul;
+ TH2* fhTrklTyTx;
TH2* fhTrklTtHMul;
TH2* fhTrklVelHMul;
TH2* fhTrklT0HMul;
@@ -247,10 +266,12 @@ private:
std::vector<TH1*> vhPullZ;
std::vector<TH1*> vhPullT;
std::vector<TH1*> vhPullTB;
+ std::vector<TH1*> vhTrefRms;
std::vector<TH2*> vhResidualTBWalk;
std::vector<TH2*> vhResidualYWalk;
std::vector<TH2*> vhXY_AllTracks; // for monitoring
std::vector<TH2*> vhXY_AllStations; // for efficiency estimation
+ std::vector<TH2*> vhXY_AllFitStations; // for efficiency estimation
std::vector<TH2*> vhXY_MissedStation; // for efficiency estimation
std::vector<TH3*> vhXY_DX;
std::vector<TH3*> vhXY_DY;
@@ -264,8 +285,9 @@ private:
TH2* fhVTX_XY0;
TH2* fhVTX_DT0_Norm;
- Int_t fTypeStation[100]; // FIXME fixed array size
- TString fOutHstFileName; // name of the histogram output file name with Calibration Parameters
+ Int_t fTypeStation[1000]; // FIXME fixed array size
+ std::vector<int> fiStationStatus; // counter status with Geo index in calibration process
+ TString fOutHstFileName; // name of the histogram output file name with Calibration Parameters
Bool_t LoadCalParameter();
Bool_t WriteHistos();
@@ -292,7 +314,8 @@ private:
Int_t fiCorMode;
Int_t fiBeamCounter;
Int_t fiStationMaxHMul;
- Double_t fTtTarg; // expected average slope in ps/cm
+ Double_t fTtTarg; // expected average slope 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
Double_t fVTX_X; // measured event wise vertex x
@@ -315,6 +338,7 @@ private:
Double_t fdRefVelMean;
Double_t fdRefDVel;
Double_t fdR0Lim;
+ Double_t fdTtMin;
TTimeStamp fStart;
TTimeStamp fStop;
@@ -323,7 +347,8 @@ private:
Double_t fdBeamMomentumLab; // beam momentum in lab frame [AGeV/c]
Bool_t fbRemoveSignalPropagationTime;
Int_t fiBeamMaxHMul;
- Int_t fiCalOpt;
+ int fiCalOpt;
+ Bool_t fbDoHistos = kTRUE;
ClassDef(CbmTofFindTracks, 1);
};
diff --git a/reco/detectors/tof/CbmTofTrackFinderNN.cxx b/reco/detectors/tof/CbmTofTrackFinderNN.cxx
index 91ce00f2179693ecba7c7c5603fa6d7f527af9fb..edaa1f80d8306ce1d08a3ace0727e85cea4b31b9 100644
--- a/reco/detectors/tof/CbmTofTrackFinderNN.cxx
+++ b/reco/detectors/tof/CbmTofTrackFinderNN.cxx
@@ -98,6 +98,8 @@ CbmTofTrackFinderNN::CbmTofTrackFinderNN(const CbmTofTrackFinderNN& finder)
, fPosYMaxScal(0.55)
, fTracks()
, fvTrkVec()
+ , fiAddVertex(0)
+ , fiVtxNbTrksMin(0)
{
// action
fHits = finder.fHits;
@@ -133,6 +135,9 @@ void CbmTofTrackFinderNN::Init()
fMinuit.Initialize();
LOG(info) << "MaxTofTimeDifference = " << fMaxTofTimeDifference;
+ if (fiAddVertex > 0) LOG(info) << "AddVertex() will be used with option " << fiAddVertex;
+ else
+ LOG(info) << "AddVertex() will not be used";
}
Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTracks)
@@ -142,25 +147,6 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
fOutTracks = fTofTracks; //new TClonesArray("CbmTofTracklet");
//fTracks = new TClonesArray("CbmTofTracklet");
//if (0 == fFindTracks->GetStationType(0)){ // Generate Pseudo TofHit at origin
- if (0) { // == fFindTracks->GetAddrOfStation(0)) { // generate new track seed, disabled
- fFindTracks->SetStation(0, 0, 0, 0);
- const Int_t iDetId = CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 0);
- const TVector3 hitPos(0., 0., 0.);
- const TVector3 hitPosErr(0.5, 0.5, 0.5); // initialize fake hit error
- const Double_t dTime0 = 0.; // FIXME
-
- Int_t iNbHits = fHits->GetEntriesFast();
- /*CbmTofHit *pHit0 =*/new ((*fHits)[iNbHits])
- CbmTofHit(iDetId, hitPos,
- hitPosErr, // local detector coordinates
- iNbHits, // this number is used as reference!!
- dTime0, // Time of hit
- 0, //vPtsRef.size(), // flag = number of TofPoints generating the cluster
- 0);
- LOG(debug1) << "CbmTofTrackFinderNN::DoFind: Fake Hit at origin added at position " << iNbHits
- << Form(", DetId 0x%08x", iDetId);
- }
-
// fvTrkMap.resize(fHits->GetEntriesFast());
fvTrkVec.resize(fHits->GetEntriesFast());
LOG(debug2) << "<I> TrkMap/Vec resized for " << fHits->GetEntriesFast() << " entries ";
@@ -169,7 +155,7 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
fvTrkVec[iHit].clear();
}
- LOG(debug) << "MaxTofTimeDifference = " << fMaxTofTimeDifference;
+ LOG(debug1) << "MaxTofTimeDifference = " << fMaxTofTimeDifference;
Int_t iNTrks = 0;
Int_t iSt0 = -1;
Int_t iSt1 = 0;
@@ -184,15 +170,15 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
// Int_t iSmType = CbmTofAddress::GetSmType( iAddr ); (VF) not used
if (HitUsed(iHit) == 1 && iAddr != fFindTracks->GetBeamCounter())
continue; // skip used Hits except for BeamCounter
- LOG(debug1) << Form("<I> TofTracklet Chkseed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
+ LOG(debug2) << Form("<I> TofTracklet Chkseed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
"0x%08x - X %6.2f, Y %6.2f Z %6.2f R %6.2f T %6.2f TM %lu",
iSt0, iSt1, fiNtrks, iHit, pHit->GetAddress(), pHit->GetX(), pHit->GetY(), pHit->GetZ(),
pHit->GetR(), pHit->GetTime(), fvTrkVec[iHit].size());
if (iAddr == fFindTracks->GetAddrOfStation(iSt0)) { // generate new track seed
- LOG(debug) << Form("<I> TofTracklet seed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
- "0x%08x - X %6.2f, Y %6.2f Z %6.2f T %6.2f TM %lu",
- iSt0, iSt1, fiNtrks, iHit, pHit->GetAddress(), pHit->GetX(), pHit->GetY(), pHit->GetZ(),
- pHit->GetTime(), fvTrkVec[iHit].size());
+ LOG(debug1) << Form("<I> TofTracklet seed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
+ "0x%08x - X %6.2f, Y %6.2f Z %6.2f T %6.2f TM %lu",
+ iSt0, iSt1, fiNtrks, iHit, pHit->GetAddress(), pHit->GetX(), pHit->GetY(), pHit->GetZ(),
+ pHit->GetTime(), fvTrkVec[iHit].size());
Int_t iChId = pHit->GetAddress();
CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
@@ -203,8 +189,7 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
if (1) { // iSmType>0) { // prevent geometry inspection for FAKE hits
if (NULL == fChannelInfo) {
- LOG(fatal) << "<D> CbmTofTrackFinderNN::DoFind0: Invalid Channel "
- "Pointer for ChId "
+ LOG(fatal) << "<D> CbmTofTrackFinderNN::DoFind0: Invalid Channel Pointer for ChId "
<< Form(" 0x%08x ", iChId) << ", Ch " << iCh;
// continue;
}
@@ -219,12 +204,12 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
gGeoManager->MasterToLocal(hitpos, hitpos_local);
}
dSizey = fChannelInfo->GetSizey();
- LOG(debug) << Form("<D> TofTracklet start %d, Hit %d - yloc %6.2f, dy %6.2f, Scal "
- "%6.2f -> stations 0x%08x, 0x%08x,",
- fiNtrks, iHit, hitpos_local[1], dSizey, fPosYMaxScal,
- fFindTracks->GetAddrOfStation(iSt0), fFindTracks->GetAddrOfStation(iSt1));
}
}
+ LOG(debug1) << Form(
+ "<D> Tracklet start %d, Hit %d - yloc %6.2f, dy %6.2f, scal %6.2f -> addrs 0x%08x, 0x%08x", fiNtrks, iHit,
+ hitpos_local[1], dSizey, fPosYMaxScal, fFindTracks->GetAddrOfStation(iSt0),
+ fFindTracks->GetAddrOfStation(iSt1));
if (TMath::Abs(hitpos_local[1]) < dSizey * fPosYMaxScal)
for (Int_t iHit1 = 0; iHit1 < fHits->GetEntriesFast(); iHit1++) // loop over all Hits (order unknown)
@@ -242,9 +227,9 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
Double_t hitpos1_local[3] = {3 * 0.};
Double_t dSizey1 = 1.;
if (NULL == fChannelInfo1) {
- LOG(debug) << "CbmTofTrackFinderNN::DoFindi: Invalid Channel "
- "Pointer for ChId "
- << Form(" 0x%08x ", iChId1) << ", Ch " << iCh1;
+ LOG(debug1) << "CbmTofTrackFinderNN::DoFindi: Invalid Channel "
+ "Pointer for ChId "
+ << Form(" 0x%08x ", iChId1) << ", Ch " << iCh1;
// continue;
}
else {
@@ -266,19 +251,21 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
Double_t dDist = TMath::Sqrt(TMath::Power((pHit->GetX() - pHit1->GetX()), 2)
+ TMath::Power((pHit->GetY() - pHit1->GetY()), 2)
+ TMath::Power((pHit->GetZ() - pHit1->GetZ()), 2));
- LOG(debug1) << Form("<I> TofTracklet %d, Hits %d, %d, add = 0x%08x,0x%08x - DT "
- "%6.2f, Tx %6.2f Ty %6.2f Tt %6.2f pos %6.2f %6.2f %6.2f ",
- fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress(), dDT, dTx, dTy,
+ LOG(debug2) << Form("<I> TofTracklet %d, Hits %d, %d, add = 0x%08x,0x%08x - DT "
+ "%6.2f, DD %6.2f, Tx %6.2f Ty %6.2f Tt %6.2f pos %6.2f %6.2f %6.2f ",
+ fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress(), dDT, dDist, dTx, dTy,
dDT / dLz, hitpos1_local[0], hitpos1_local[1], hitpos1_local[2]);
- LOG(debug3) << Form(" selection: y %6.2f < %6.2f, T %6.2f < %6.2f, "
- "dTpos %6.2f < %6.2f, Abs(%6.2f - %6.2f) < %6.2f",
+ LOG(debug3) << Form(" DT %f from %f - %f ", dDT, pHit1->GetTime(), pHit->GetTime());
+ LOG(debug3) << Form(" selection: y %6.2f < %6.2f, T %f < %6.5f, "
+ "Tx Abs(%6.2f - %6.2f) < %6.2f, Ty Abs(%6.2f - %6.2f) < %6.2f",
hitpos1_local[1], dSizey1 * fPosYMaxScal, dDT / dLz, fMaxTofTimeDifference, dTx,
- fTxLIM, dTy, fTyMean, fTyLIM);
+ fTxMean, fTxLIM, dTy, fTyMean, fTyLIM);
if (TMath::Abs(hitpos1_local[1]) < dSizey1 * fPosYMaxScal)
if (TMath::Abs(dDT / dLz) < fMaxTofTimeDifference
- && TMath::Abs(dDT / dDist) > 0.005 // FIXME: numeric constant in code
+ //&& TMath::Abs(dDT / dDist) > 0.0001 // FIXME: numeric constant in code
&& TMath::Abs(dTx - fTxMean) < fTxLIM && TMath::Abs(dTy - fTyMean) < fTyLIM) {
+ LOG(debug3) << "Construct new tracklet";
CbmTofTracklet* pTrk = new CbmTofTracklet();
fTracks.push_back(pTrk);
pTrk->SetTofHitIndex(iHit, iAddr, pHit); // store 1. Hit index
@@ -313,11 +300,11 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
tPar->SetTx(dTx);
tPar->SetTy(dTy);
- LOG(debug) << Form("<I> TofTracklet %d, Hits %d, %d initialized, "
- "add 0x%08x,0x%08x, DT %6.3f, Sgn %2.0f, DR "
- "%6.3f, T0 %6.2f, Tt %6.4f ",
- fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress(), dDT, dSign, dR,
- pTrk->GetT0(), pTrk->GetTt())
+ LOG(debug1) << Form("<I> TofTracklet %d, Hits %d, %d initialized, "
+ "add 0x%08x,0x%08x, DT %6.3f, Sgn %2.0f, DR "
+ "%6.3f, T0 %6.2f, Tt %6.4f ",
+ fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress(), dDT, dSign, dR,
+ pTrk->GetT0(), pTrk->GetTt())
// << tPar->ToString()
;
}
@@ -358,9 +345,9 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
Double_t dSizey = 1.;
if (NULL == fChannelInfo) {
- LOG(debug) << "CbmTofTrackFinderNN::DoFind: Invalid Channel "
- "Pointer from Hit "
- << iHit << " for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
+ LOG(debug1) << "CbmTofTrackFinderNN::DoFind: Invalid Channel "
+ "Pointer from Hit "
+ << iHit << " for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
// continue;
}
else {
@@ -394,33 +381,37 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
// pTrk->GetFitT(pHit->GetR());
Double_t dChi =
- TMath::Sqrt((TMath::Power(TMath::Abs(dTex - pHit->GetTime()) / fFindTracks->GetSigT(iAddr), 2)
- + TMath::Power(TMath::Abs(dXex - pHit->GetX()) / fFindTracks->GetSigX(iAddr), 2)
- + TMath::Power(TMath::Abs(dYex - pHit->GetY()) / fFindTracks->GetSigY(iAddr), 2))
+ TMath::Sqrt((TMath::Power(TMath::Abs(dTex - pHit->GetTime()) / fFindTracks->GetStationSigT(iDet), 2)
+ + TMath::Power(TMath::Abs(dXex - pHit->GetX()) / fFindTracks->GetStationSigX(iDet), 2)
+ + TMath::Power(TMath::Abs(dYex - pHit->GetY()) / fFindTracks->GetStationSigY(iDet), 2))
/ 3);
- LOG(debug1) << Form("<IP> TofTracklet %lu, HMul %d, Hits %d, %d check %d, Station "
- "0x%08x: DT %f, DX %f, DY %f, Chi %f",
- iTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, iAddr,
- (dTex - pHit->GetTime()) / fFindTracks->GetSigT(iAddr),
- (dXex - pHit->GetX()) / fFindTracks->GetSigX(iAddr),
- (dYex - pHit->GetY()) / fFindTracks->GetSigY(iAddr), dChi);
+ LOG(debug2) << Form("<IP> TofTracklet %lu, HMul %d, Hits %d, %d check %d, Station "
+ "%d: DT %f, DX %f, DY %f, Chi %f",
+ iTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, iDet,
+ (dTex - pHit->GetTime()) / fFindTracks->GetStationSigT(iDet),
+ (dXex - pHit->GetX()) / fFindTracks->GetStationSigX(iDet),
+ (dYex - pHit->GetY()) / fFindTracks->GetStationSigY(iDet), dChi);
+
+ if (fFindTracks->GetStationSigT(iDet) == 0 || fFindTracks->GetStationSigX(iDet) == 0
+ || fFindTracks->GetStationSigY(iDet) == 0)
+ LOG(fatal) << Form("Missing resolution for station %d, addr 0x%08x", iDet, iAddr);
if (dChi < fSIGLIM // FIXME: should scale limit with material budget between hit and track reference
* (pTrk->GetNofHits() < fFindTracks->GetNReqStations() - 1
? 1.
: fSIGLIMMOD)) { // extend and update tracklet
- LOG(debug) << Form("<IP> TofTracklet %lu, HMul %d, Hits %d, %d "
- "mark for extension by %d, add = 0x%08x, DT "
- "%6.2f, DX %6.2f, DY=%6.2f ",
- iTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, pHit->GetAddress(),
- dTex - pHit->GetTime(), dXex - pHit->GetX(), dYex - pHit->GetY())
- << tPar->ToString();
+ LOG(debug1) << Form("<IP> TofTracklet %lu, HMul %d, Hits %d, %d "
+ "mark for extension by %d, add = 0x%08x, DT "
+ "%6.2f, DX %6.2f, DY=%6.2f ",
+ iTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, pHit->GetAddress(),
+ dTex - pHit->GetTime(), dXex - pHit->GetX(), dYex - pHit->GetY())
+ << tPar->ToString();
if (iNCand > 0) {
- LOG(debug) << Form("CbmTofTrackFinderNN::DoFind new match %d "
- "of Hit %d, Trk %lu, chi2 = %f",
- iNCand, iHit, iTrk, dChi);
+ LOG(debug1) << Form("CbmTofTrackFinderNN::DoFind new match %d "
+ "of Hit %d, Trk %lu, chi2 = %f",
+ iNCand, iHit, iTrk, dChi);
iNCand++;
if (iNCand == MAXNCAND) iNCand--; // Limit index to maximum
@@ -435,15 +426,15 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
iHitInd[iCand] = iHit;
dChi2[iCand] = dChi;
dChi2[iNCand] = 1.E8;
- LOG(debug1) << Form(" <D> candidate inserted at pos %d", iCand);
+ LOG(debug2) << Form(" <D> candidate inserted at pos %d", iCand);
break;
}
}
}
else {
- LOG(debug) << Form("CbmTofTrackFinderNN::DoFind first match "
- "%d of Hit %d, Trk %p, chi2 = %f",
- iNCand, iHit, pTrk, dChi);
+ LOG(debug1) << Form("CbmTofTrackFinderNN::DoFind first match "
+ "%d of Hit %d, Trk %p, chi2 = %f",
+ iNCand, iHit, pTrk, dChi);
pTrkInd[iNCand] = pTrk;
iHitInd[iNCand] = iHit;
dChi2[iNCand] = dChi; // relative quality measure
@@ -459,14 +450,15 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
while (iNCand > 0) { // at least one matching hit - trk pair found
CbmTofTracklet* pTrk = pTrkInd[0];
if (NULL == pTrk) continue;
- LOG(debug) << Form("%d hit match candidates in station %d to %lu TofTracklets", iNCand, iDet, fTracks.size());
+ LOG(debug1) << Form("%d hit match candidates in station %d to %lu TofTracklets", iNCand, iDet,
+ fTracks.size());
for (Int_t iM = 0; iM < iNCand; iM++) {
pTrk = (CbmTofTracklet*) pTrkInd[iM];
if (NULL == pTrk) break;
std::vector<CbmTofTracklet*>::iterator it = std::find(fTracks.begin(), fTracks.end(), pTrk);
if (it == fTracks.end()) break; // track candidate not existing
- LOG(debug1) << "\t"
+ LOG(debug2) << "\t"
<< Form("Hit %d, Trk %p with chi2 %f (%f)", iHitInd[iM], pTrkInd[iM], dChi2[iM],
pTrk->GetMatChi2(fFindTracks->GetAddrOfStation(iDet)));
}
@@ -477,7 +469,7 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
size_t iTr = 0;
for (; iTr < fTracks.size(); iTr++) {
if (fTracks[iTr] == pTrk) {
- LOG(debug) << "Track " << pTrk << " active at pos " << iTr;
+ LOG(debug1) << "Track " << pTrk << " active at pos " << iTr;
break;
}
}
@@ -501,7 +493,7 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
CbmTofHit* pHit = (CbmTofHit*) fHits->At(iHit);
Int_t iAddr = (pHit->GetAddress() & DetMask);
if (Double_t dLastChi2 = pTrk->GetMatChi2(fFindTracks->GetAddrOfStation(iDet)) == -1.) {
- LOG(debug1) << Form(" -D- Add hit %d at %p, Addr 0x%08x, Chi2 %6.2f to size %u", iHit, pHit, iAddr,
+ LOG(debug2) << Form(" -D- Add hit %d at %p, Addr 0x%08x, Chi2 %6.2f to size %u", iHit, pHit, iAddr,
dChi2[0], pTrk->GetNofHits());
pTrk->AddTofHitIndex(iHit, iAddr, pHit,
dChi2[0]); // store next Hit index with matching chi2
@@ -509,10 +501,11 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
fvTrkVec[iHit].push_back(pTrk);
Line3Dfit(pTrk); // full MINUIT fit overwrites ParamLast!
Bool_t bkeep = kFALSE;
- if (pTrk->GetChiSq() > fChiMaxAccept) {
- LOG(debug) << Form("Add hit %d invalidates tracklet with Chi "
- "%6.2f > %6.2f -> undo ",
- iHit, pTrk->GetChiSq(), fChiMaxAccept);
+ if ((pTrk->GetChiSq() > fChiMaxAccept)
+ || (fFindTracks->GetR0Lim() > 0 && pTrk->GetR0() > fFindTracks->GetR0Lim())) {
+ LOG(debug1) << Form("Add hit %d invalidates tracklet with Chi "
+ "%6.2f > %6.2f -> undo ",
+ iHit, pTrk->GetChiSq(), fChiMaxAccept);
fvTrkVec[iHit].pop_back();
pTrk->RemoveTofHitIndex(iHit, iAddr, pHit, dChi2[0]);
Line3Dfit(pTrk); //restore old status
@@ -550,8 +543,8 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
// TODO remove tracklet assigment of old hit! FIXME
}
else {
- LOG(debug) << Form(" -D- Ignore %d, Det %d, Addr 0x%08x, at 0x%p, Chi2 %6.2f", iHit, iDet, iAddr, pHit,
- dChi2[0]);
+ LOG(debug1) << Form(" -D- Ignore %d, Det %d, Addr 0x%08x, at 0x%p, Chi2 %6.2f", iHit, iDet, iAddr,
+ pHit, dChi2[0]);
// Form new seeds
//if (iDet<(fFindTracks->GetNStations()-1)) TrklSeed(fHits,fTracks,iHit);
break;
@@ -580,15 +573,15 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
// update inverse velocity
Double_t dTt = pTrk->GetTt();
- LOG(debug) << Form("<Res> TofTracklet %p, HMul %d, Hits %d, %d, %d, "
- "NDF %d, Chi2 %6.2f, T0 %6.2f, Tt %6.4f ",
- pTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, pTrk->GetNDF(), pTrk->GetChiSq(),
- pTrk->GetTime(), dTt);
+ LOG(debug1) << Form("<Res> TofTracklet %p, HMul %d, Hits %d, %d, %d, "
+ "NDF %d, Chi2 %6.2f, T0 %6.2f, Tt %6.4f ",
+ pTrk, pTrk->GetNofHits(), iHit0, iHit1, iHit, pTrk->GetNDF(), pTrk->GetChiSq(),
+ pTrk->GetTime(), dTt);
PrintStatus((char*) "<Res> ");
- LOG(debug1) << " Match loop status: NCand " << iNCand << ", iDet " << iDet;
+ LOG(debug2) << " Match loop status: NCand " << iNCand << ", iDet " << iDet;
- /* live display insert
- if(fair::Logger::Logging(fair::Severity::debug3)) // update event display, if initialized
+ /* live display insert
+ if(fair::Logger::Logging(fair::Severity::debug3)) // update event display, if initialized
{
Int_t ii;
CbmEvDisTracks* fDis = CbmEvDisTracks::Instance();
@@ -604,7 +597,7 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
}
cout << " fDis "<<fDis<<" with "<<fiNtrks<<" tracks, to continue type 0 ! "<<endl;
scanf("%d",&ii);
- } // end of live display
+ } // end of live display
*/
} // end of while(iNCand>0)
} // detector loop (propagate) end
@@ -613,9 +606,15 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
//fTracks->Compress();
//fTofTracks = fTracks;
- // copy fTracks -> fTofTracks / fOutTracks
+ //fFindTracks->PrintSetup();
+ // Add Vertex as additional point
+ //if(fbAddVertex) LOG(fatal)<<"Debugging step found AddVertex";
+ if (fiAddVertex > 0) AddVertex();
+ //fFindTracks->PrintSetup();
+ // copy fTracks -> fTofTracks / fOutTracks
+ LOG(debug1) << "Clean-up " << fTracks.size() << " tracklet candidates";
fiNtrks = 0;
for (size_t iTr = 0; iTr < fTracks.size(); iTr++) {
if (fTracks[iTr] == NULL) continue;
@@ -624,22 +623,25 @@ Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits, TClonesArray* fTofTrac
CbmTofTracklet* pTrk = new ((*fTofTracks)[fiNtrks++]) CbmTofTracklet(*fTracks[iTr]);
if (fair::Logger::Logging(fair::Severity::debug)) {
- LOG(info) << "Found Trkl " << iTr << ", ";
+ LOG(info) << "Found Trkl " << iTr;
pTrk->PrintInfo();
}
for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) { // mark used Hit
CbmTofHit* pHit = (CbmTofHit*) fHits->At(pTrk->GetHitIndex(iHit));
pHit->SetFlag(pHit->GetFlag() + 100.);
- LOG(debug) << Form(" hit %d at %d flagged to %d ", iHit, pTrk->GetHitIndex(iHit), pHit->GetFlag());
+ LOG(debug1) << Form(" hit %d at %d flagged to %d ", iHit, pTrk->GetHitIndex(iHit), pHit->GetFlag());
}
}
- PrintStatus((char*) "<D> Final result");
+ if (fair::Logger::Logging(fair::Severity::debug)) {
+ fFindTracks->PrintSetup();
+ PrintStatus((char*) "<D> Final result");
+ }
for (size_t iTr = 0; iTr < fTracks.size(); iTr++) {
if (fTracks[iTr] == NULL) continue;
fTracks[iTr]->Delete();
//delete fTracks[iTr];
- LOG(debug) << Form("<I> TofTracklet %lu, %p deleted", iTr, fTracks[iTr]);
+ LOG(debug1) << Form("<I> TofTracklet %lu, %p deleted", iTr, fTracks[iTr]);
}
fTracks.resize(0); //cleanup
// fFindTracks->PrintSetup();
@@ -670,8 +672,8 @@ void CbmTofTrackFinderNN::TrklSeed(Int_t iHit)
Double_t dSizey1 = 1.;
if (NULL == fChannelInfo1) {
- LOG(debug) << "CbmTofTrackFinderNN::DoFindp: Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId1)
- << ", Ch " << iCh1;
+ LOG(debug1) << "CbmTofTrackFinderNN::DoFindp: Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId1)
+ << ", Ch " << iCh1;
// continue;
}
else {
@@ -689,7 +691,7 @@ void CbmTofTrackFinderNN::TrklSeed(Int_t iHit)
Double_t dTx = (pHit->GetX() - pHit1->GetX()) / dLz;
Double_t dTy = (pHit->GetY() - pHit1->GetY()) / dLz;
Int_t iUsed = HitUsed(iHit1);
- LOG(debug1) << Form("<ISeed> TofTracklet %d, Hits %d, %d, used %d check, add = "
+ LOG(debug2) << Form("<ISeed> TofTracklet %d, Hits %d, %d, used %d check, add = "
"0x%08x,0x%08x - DT %6.2f, Tx %6.2f Ty %6.2f ",
fiNtrks, iHit, iHit1, iUsed, pHit->GetAddress(), pHit1->GetAddress(), dDT, dTx, dTy);
if (TMath::Abs(hitpos1_local[1]) < dSizey1 * fPosYMaxScal)
@@ -724,9 +726,9 @@ void CbmTofTrackFinderNN::TrklSeed(Int_t iHit)
tPar->SetLz(dLz);
tPar->SetTx(dTx);
tPar->SetTy(dTy);
- LOG(debug) << Form("<DSeed> TofTracklet %d, Hits %d, %d add "
- "initialized, add = 0x%08x,0x%08x ",
- fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress());
+ LOG(debug1) << Form("<DSeed> TofTracklet %d, Hits %d, %d add "
+ "initialized, add = 0x%08x,0x%08x ",
+ fiNtrks, iHit, iHit1, pHit->GetAddress(), pHit1->GetAddress());
PrintStatus((char*) "after DSeed");
}
}
@@ -783,12 +785,19 @@ void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk)
<< Form(", %d.Hit(%d) at ind %d with %d(%d) registered tracks", iHit, pTrk->GetNofHits(), iHitInd,
(int) fvTrkVec[iHitInd].size(), NTrks);
//if(fvTrkVec[iHitInd].size()==1) break;
+ Int_t iTrkPos = 0;
for (std::vector<CbmTofTracklet*>::iterator iT = fvTrkVec[iHitInd].begin(); iT != fvTrkVec[iHitInd].end();
iT++) {
iterClean = 0;
+ LOG(debug2) << "Inspect track " << iTrkPos << " of " << fvTrkVec[iHitInd].size() << " for HitInd " << iHitInd;
+ if (iTrkPos >= (Int_t) fvTrkVec[iHitInd].size()) break;
+ iTrkPos++;
if (!Active(*iT)) break; // check whether tracklet is still active
LOG(debug2) << " <D2> process Trk " << *iT << " with " << (*iT)->GetNofHits() << " hits";
-
+ if (*iT == pTrk) {
+ LOG(debug2) << " <D2a> continue with next tracklet ";
+ continue;
+ }
for (Int_t iH = 0; iH < (*iT)->GetNofHits(); iH++) {
if (!Active(*iT)) break; // check whether tracklet is still active
Int_t iHi = (*iT)->GetTofHitIndex(iH);
@@ -815,7 +824,6 @@ void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk)
<< iterClean << ", hit " << iHi << ", size " << fvTrkVec[iHi].size();
if (*it != pTrk) {
size_t iTr = 0;
- ;
for (iTr = 0; iTr < fTracks.size(); iTr++) {
if (fTracks[iTr] == *it) {
LOG(debug2) << Form(" found track entry %p(%d) at %lu "
@@ -874,24 +882,24 @@ void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk)
LOG(debug2) << "<D8> remove tracklet at pos " << iTr;
fTracks.erase(fTracks.begin() + iTr);
fiNtrks--;
-
- LOG(debug2) << "Erase1 for pTrk " << pTrk << ", at " << iTr << ", hit " << iHi << ", size "
- << fvTrkVec[iHi].size();
-
- PrintStatus((char*) "UpdateTrackList::Erase1");
+ if (fair::Logger::Logging(fair::Severity::debug2)) {
+ LOG(debug2) << "Erase1 for pTrk " << pTrk << ", at " << iTr << ", hit " << iHi << ", size "
+ << fvTrkVec[iHi].size();
+ PrintStatus((char*) "UpdateTrackList::Erase1");
+ }
/*
if(fvTrkVec[iHi].size() == 1) {
fvTrkVec[iHi].clear();
LOG(debug2) << " clear1 for pTrk "<<pTrk<<", hit "<<iHi<<", size "<<fvTrkVec[iHi].size()
;
- goto loopclean;
+ goto loopclean;
}else{
it=fvTrkVec[iHi].erase(it);
//NTrks--;
LOG(debug2) << " erase3 for "<<iTrk<<" at "<<pTrk<<", hit "<<iHi
<<", size "<<fvTrkVec[iHi].size()<<", "<<NTrks
- ;
+ ;
}
*/
@@ -905,9 +913,13 @@ void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk)
// if(pTrk == *iT) goto loopclean; //
}
} // end of loop over tracks referenced by hit iHi
+ LOG(debug2) << "Track loop of iHi = " << iHi << " finished";
}
+ if (!Active(*iT)) break; // check whether tracklet is still active
}
+ LOG(debug2) << "Hit loop of iTrkPos = " << iTrkPos << " finished";
}
+ LOG(debug2) << "Track loop of iHitInd = " << iHitInd << " finished";
}
}
}
@@ -915,8 +927,8 @@ void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk)
void CbmTofTrackFinderNN::PrintStatus(char* cComment)
{
- LOG(debug) << Form("<PS %s> for fiNtrks = %d tracks out of %d fTracks.size() ", cComment, fiNtrks,
- (int) fTracks.size());
+ LOG(debug1) << Form("<PS %s> for fiNtrks = %d tracks of %d, %d assigned hits ", cComment, fiNtrks,
+ (int) fTracks.size(), (int) fvTrkVec.size());
for (size_t it = 0; it < fTracks.size(); it++) {
CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[it];
@@ -930,11 +942,12 @@ void CbmTofTrackFinderNN::PrintStatus(char* cComment)
}
sTrk += Form(", ChiSq %7.1f", pTrk->GetChiSq());
sTrk += Form(", Tt %6.4f", pTrk->GetTt());
- LOG(debug) << sTrk;
+ LOG(debug1) << sTrk;
}
for (size_t ih = 0; ih < fvTrkVec.size(); ih++) {
CbmTofHit* pHit = (CbmTofHit*) fHits->At(ih);
+ if (NULL == pHit) LOG(fatal) << "<E> missing pointer for hit " << ih;
Int_t iAddr = (pHit->GetAddress() & DetMask);
Int_t iSt = fFindTracks->GetStationOfAddr(iAddr);
TString sTrk = "";
@@ -945,8 +958,8 @@ void CbmTofTrackFinderNN::PrintStatus(char* cComment)
CbmTofTracklet* pTrk = fvTrkVec[ih][it];
sTrk += Form(" %p ", pTrk);
}
- LOG(debug) << sTrk;
}
+ LOG(debug1) << sTrk;
}
}
@@ -961,12 +974,21 @@ Bool_t CbmTofTrackFinderNN::Active(CbmTofTracklet* pCheck)
}
void CbmTofTrackFinderNN::Line3Dfit(CbmTofTracklet* pTrk)
+{
+ Int_t iDetAddr = 0;
+ Line3Dfit(pTrk, iDetAddr);
+}
+
+void CbmTofTrackFinderNN::Line3Dfit(CbmTofTracklet* pTrk, Int_t iDetAddr)
{
TGraph2DErrors* gr = new TGraph2DErrors();
// Fill the 2D graph
// generate graph with the 3d points
+ Int_t Np = 0;
for (Int_t N = 0; N < pTrk->GetNofHits(); N++) {
+ if (((pTrk->GetTofHitPointer(N))->GetAddress() & DetMask) == iDetAddr) continue; //skip specific detector
+ Np++;
double x, y, z = 0;
x = (pTrk->GetTofHitPointer(N))->GetX();
y = (pTrk->GetTofHitPointer(N))->GetY();
@@ -977,8 +999,8 @@ void CbmTofTrackFinderNN::Line3Dfit(CbmTofTracklet* pTrk)
dy = (pTrk->GetTofHitPointer(N))->GetDy();
dz = (pTrk->GetTofHitPointer(N))->GetDz(); //FIXME
gr->SetPointError(N, dx, dy, dz);
- LOG(debug) << "Line3Dfit add N = " << N << ",\t" << pTrk->GetTofHitIndex(N) << ",\t" << x << ",\t" << y << ",\t"
- << z << ",\t" << dx << ",\t" << dy << ",\t" << dz;
+ LOG(debug2) << "Line3Dfit add N = " << N << ",\t" << pTrk->GetTofHitIndex(N) << ",\t" << x << ",\t" << y << ",\t"
+ << z << ",\t" << dx << ",\t" << dy << ",\t" << dz;
}
// fit the graph now
Double_t pStart[4] = {0., 0., 0., 0.};
@@ -986,25 +1008,137 @@ void CbmTofTrackFinderNN::Line3Dfit(CbmTofTracklet* pTrk)
pStart[1] = (pTrk->GetTrackParameter())->GetTx();
pStart[2] = pTrk->GetFitY(0.);
pStart[3] = (pTrk->GetTrackParameter())->GetTy();
- LOG(debug) << "Line3Dfit init: X0 " << pStart[0] << ", TX " << pStart[1] << ", Y0 " << pStart[2] << ", TY "
- << pStart[3];
+ LOG(debug2) << "Line3Dfit init: X0 " << pStart[0] << ", TX " << pStart[1] << ", Y0 " << pStart[2] << ", TY "
+ << pStart[3];
fMinuit.DoFit(gr, pStart);
//gr->Draw("err p0");
gr->Delete();
Double_t* dRes;
dRes = fMinuit.GetParFit();
- LOG(debug) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0] << ", TX " << dRes[1] << ", Y0 "
- << dRes[2] << ", TY " << dRes[3] << ", Chi2DoF: " << fMinuit.GetChi2DoF();
+ LOG(debug2) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0] << ", TX " << dRes[1] << ", Y0 "
+ << dRes[2] << ", TY " << dRes[3] << ", Chi2DoF: " << fMinuit.GetChi2DoF();
(pTrk->GetTrackParameter())->SetX(dRes[0]);
(pTrk->GetTrackParameter())->SetY(dRes[2]);
(pTrk->GetTrackParameter())->SetZ(0.);
(pTrk->GetTrackParameter())->SetTx(dRes[1]);
(pTrk->GetTrackParameter())->SetTy(dRes[3]);
- (pTrk->GetTrackParameter())->SetQp(1.); // FIXME
- pTrk->SetChiSq(fMinuit.GetChi2DoF() / pTrk->GetNofHits());
+ (pTrk->GetTrackParameter())->SetQp(1.); // FIXME
+ pTrk->SetChiSq(fMinuit.GetChi2DoF() / Np); //pTrk->GetNofHits());
// empirical to equilibrate bias on hit multiplicity!!!
}
+void CbmTofTrackFinderNN::AddVertex()
+{
+ if (fTracks.size() < 2) return;
+
+ Double_t dTvtx = 0.;
+ Double_t dT2vtx = 0.;
+ Double_t dSigTvtx = 0.;
+ Double_t nValidTracks = 0.;
+ Double_t dToff = 0.;
+ double dVx = 0.;
+ double dVy = 0.;
+ double dVx2 = 0.;
+ double dVy2 = 0.;
+ for (size_t it = 0; it < fTracks.size(); it++) {
+ CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[it];
+ if (NULL == pTrk) continue;
+ if (pTrk->GetR0() < fFindTracks->GetR0Lim()) {
+ if (nValidTracks == 0) dToff = pTrk->GetT0();
+ dTvtx += pTrk->GetT0() - dToff;
+ dVx += pTrk->GetFitX(0.);
+ dVy += pTrk->GetFitY(0.);
+ dVx2 += dVx * dVx;
+ dVy2 += dVy * dVy;
+ dT2vtx += (pTrk->GetT0() - dToff) * (pTrk->GetT0() - dToff);
+ nValidTracks += 1;
+ }
+ }
+ LOG(debug1) << Form("AddVertex valid tracks %3.0f: %6.3f, %6.3f, %15.3f", nValidTracks, dTvtx, dT2vtx, dToff);
+ if (nValidTracks >= fiVtxNbTrksMin) { // generate virtual hit
+ dTvtx /= nValidTracks;
+ dT2vtx /= nValidTracks;
+ dSigTvtx = TMath::Sqrt(dT2vtx - dTvtx * dTvtx);
+ dVx /= nValidTracks;
+ dVy /= nValidTracks;
+ dVx2 /= nValidTracks;
+ dVy2 /= nValidTracks;
+ LOG(debug1) << Form("AddVertex time %15.3f, sig %8.3f from %3.0f tracks ", dTvtx, dSigTvtx, nValidTracks);
+ if (dSigTvtx < 1. * fChiMaxAccept) // FIXME constant in code
+ {
+ const Int_t iDetId = CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 10);
+ if (fiAddVertex == 10) {
+ dVx = 0.;
+ dVy = 0.;
+ }
+ double dSigVx = TMath::Sqrt(dVx2 - dVx * dVx);
+ double dSigVy = TMath::Sqrt(dVy2 - dVy * dVy);
+ if (fiAddVertex > 1) {
+ dSigVx /= nValidTracks;
+ dSigVy /= nValidTracks;
+ }
+ const TVector3 hitPos(dVx, dVy, 0.);
+ const TVector3 hitPosErr(dSigVx, dSigVy, 0.1); // initialize fake hit error
+ const Double_t dTime0 = dTvtx + dToff; // FIXME
+ /*
+ if( fFindTracks->GetStationOfAddr(iDetId) < 0) {
+ fFindTracks->SetStation(fFindTracks->GetNStations(), 10, 0, 0);
+ fFindTracks->SetNStations(fFindTracks->GetNStations()+1);
+ LOG(warn)<<"AddVertex increased NStations to "<<fFindTracks->GetNStations();
+ }
+ */
+ Int_t iHitVtx = fHits->GetEntriesFast();
+ CbmTofHit* pHitVtx =
+ new ((*fHits)[iHitVtx]) CbmTofHit(iDetId, hitPos,
+ hitPosErr, // local detector coordinates
+ iHitVtx, // this number is used as reference!!
+ dTime0, // Time of hit
+ 0, //vPtsRef.size(), // flag = number of TofPoints generating the cluster
+ 0);
+ LOG(debug1) << "CbmTofTrackFinderNN::DoFind: Fake Vtx Hit at pos " << iHitVtx << Form(", T0 %f ", dTime0)
+ << Form(", DetId 0x%08x TSR %d%d%d ", iDetId, CbmTofAddress::GetSmType(iDetId),
+ CbmTofAddress::GetSmId(iDetId), CbmTofAddress::GetRpcId(iDetId));
+ pHitVtx->SetTimeError(dSigTvtx);
+
+ PrintStatus((char*) "InAddVertex");
+
+ fvTrkVec.resize(fHits->GetEntriesFast());
+ LOG(debug2) << "<I> TrkMap/Vec resized for " << fHits->GetEntriesFast() << " entries ";
+
+ for (size_t it = 0; it < fTracks.size(); it++) {
+ CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[it];
+ if (NULL == pTrk) continue;
+ if (pTrk->GetR0() < fFindTracks->GetR0Lim()) {
+ pTrk->AddTofHitIndex(iHitVtx, iDetId, pHitVtx);
+ fvTrkVec[iHitVtx].push_back(pTrk);
+ Line3Dfit(pTrk); // full MINUIT fit overwrites ParamLast!
+ if (pTrk->GetChiSq() > fChiMaxAccept) {
+ LOG(debug1) << Form("Add hit %d invalidates tracklet with Chi "
+ "%6.2f > %6.2f -> undo ",
+ iHitVtx, pTrk->GetChiSq(), fChiMaxAccept);
+ fvTrkVec[iHitVtx].pop_back();
+ pTrk->RemoveTofHitIndex(iHitVtx, iDetId, pHitVtx, 0.);
+ Line3Dfit(pTrk); //restore old status
+ }
+ else {
+ // update times fit
+ Double_t dTt = pTrk->GetTt();
+ LOG(debug1) << "Vtx added to track " << it << "; chi2 " << pTrk->GetChiSq() << ", Tt " << dTt;
+ }
+ }
+ } //loop on tracks finished
+ PrintStatus((char*) "PostAddVertex");
+ }
+ }
+ /*
+ // fvTrkMap.resize(fHits->GetEntriesFast());
+ // for (Int_t iHit=0; iHit<fHits->GetEntriesFast(); iHit++) { fvTrkMap[iHit].clear();}
+ for (Int_t iHit = 0; iHit < fHits->GetEntriesFast(); iHit++) {
+ fvTrkVec[iHit].clear();
+ }
+*/
+}
+
ClassImp(CbmTofTrackFinderNN)
diff --git a/reco/detectors/tof/CbmTofTrackFinderNN.h b/reco/detectors/tof/CbmTofTrackFinderNN.h
index 33c36bb40d4ad8194b322c34ec51f3f8b1388f16..514de7abcd629aa1236e98d85915481ac359d105 100644
--- a/reco/detectors/tof/CbmTofTrackFinderNN.h
+++ b/reco/detectors/tof/CbmTofTrackFinderNN.h
@@ -3,7 +3,7 @@
Authors: Norbert Herrmann [committer], Pierre-Alain Loizeau */
/**
-nh, adapt from
+nh, adapt from
* \file CbmTrdTrackFinderIdeal.h
**/
@@ -44,8 +44,8 @@ public:
void TrklSeed(Int_t iHit);
Int_t HitUsed(Int_t iHit);
- /*
- void RemoveMultipleAssignedHits(
+ /*
+ void RemoveMultipleAssignedHits(
TClonesArray* fTofHits,
Int_t iDet
);
@@ -74,9 +74,13 @@ public:
inline Double_t GetChiMaxAccept() { return fChiMaxAccept; }
static void Line3Dfit(CbmTofTracklet* pTrk);
+ static void Line3Dfit(CbmTofTracklet* pTrk, Int_t iAddr);
Bool_t Active(CbmTofTracklet* pTrk);
void PrintStatus(char* cComm);
+ void AddVertex();
+ inline void SetAddVertex(int ival) { fiAddVertex = ival; }
+ inline void SetVtxNbTrksMin(int ival) { fiVtxNbTrksMin = ival; }
//Copy constructor
CbmTofTrackFinderNN(const CbmTofTrackFinderNN& finder);
@@ -105,7 +109,8 @@ private:
std::vector<CbmTofTracklet*> fTracks; // Tracklets to which hit is assigned
//std::vector<std::map <CbmTofTracklet *, Int_t> > fvTrkMap; // Tracklets to which hit is assigned
std::vector<std::vector<CbmTofTracklet*>> fvTrkVec; // Tracklets to which hit is assigned
-
+ int fiAddVertex;
+ int fiVtxNbTrksMin;
ClassDef(CbmTofTrackFinderNN, 1);
};
diff --git a/reco/detectors/tof/CbmTofTrackletTools.cxx b/reco/detectors/tof/CbmTofTrackletTools.cxx
index 41119cf73935f8ebb00cdc3442f035057eb2fa09..43978e982214eeabb9a22bdae821b4d5419410b2 100644
--- a/reco/detectors/tof/CbmTofTrackletTools.cxx
+++ b/reco/detectors/tof/CbmTofTrackletTools.cxx
@@ -52,27 +52,28 @@ Double_t* CbmTofTrackletTools::Line3DFit(CbmTofTracklet* pTrk, Int_t iDetId)
pStart[1] = (pTrk->GetTrackParameter())->GetTx();
pStart[2] = pTrk->GetFitY(0.);
pStart[3] = (pTrk->GetTrackParameter())->GetTy();
- LOG(debug) << "Line3DFit init: X0 " << pStart[0] << ", TX " << pStart[1] << ", Y0 " << pStart[2] << ", TY "
- << pStart[3];
+ LOG(debug1) << "Line3DFit init: X0 " << pStart[0] << ", TX " << pStart[1] << ", Y0 " << pStart[2] << ", TY "
+ << pStart[3];
fMinuit.DoFit(gr, pStart);
gr->Delete();
Double_t* dRes;
// LOG(info) << "Get fit results ";
dRes = fMinuit.GetParFit();
- LOG(debug) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0] << ", TX " << dRes[1] << ", Y0 "
- << dRes[2] << ", TY " << dRes[3] << ", Chi2DoF: " << fMinuit.GetChi2DoF();
+ LOG(debug1) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0] << ", TX " << dRes[1] << ", Y0 "
+ << dRes[2] << ", TY " << dRes[3] << ", Chi2DoF: " << fMinuit.GetChi2DoF();
return dRes;
}
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.;
//LOG(info) << "FitTt: " << pTrk->GetNofHits() << " hits, exclude " << Form("0x%08x",iDetId);
- Double_t aR[pTrk->GetNofHits() - 1];
- Double_t at[pTrk->GetNofHits() - 1];
- Double_t ae[pTrk->GetNofHits() - 1];
+ Double_t aR[pTrk->GetNofHits()];
+ Double_t at[pTrk->GetNofHits()];
+ Double_t ae[pTrk->GetNofHits()];
for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
if (iDetId == pTrk->GetTofDetIndex(iHit)) continue;
if (nValidHits == 0) {
@@ -111,14 +112,20 @@ Double_t CbmTofTrackletTools::FitTt(CbmTofTracklet* pTrk, Int_t iDetId)
}
Double_t det_cov_mat =
esum * RRsum
- - Rsum * Rsum; // Determinant of inverted Covariance Matrix -> 1/det is common Faktor of Cavariance Matrix
- //fT0=(RRsum*tsum-Rsum*Rtsum)/det_cov_mat; // Best Guess for time at origin
- return (-Rsum * tsum + esum * Rtsum) / det_cov_mat; // Best guess for inverted velocity
- //fT0Err=TMath::Sqrt(RRsum/det_cov_mat); // sqrt of (0,0) in Covariance matrix -> error on fT0
- //fTtErr=TMath::Sqrt(esum/det_cov_mat); // sqrt of (1,1) in Covariance Matrix -> error on fTt
- //fT0TtCov=-Rsum/det_cov_mat; // (0,1)=(1,0) in Covariance Matrix -> cov(fT0,fTt)
-
- //fT0+=yoffset; // Adding yoffset again
+ - Rsum * Rsum; // Determinant of inverted Covariance Matrix -> 1/det is common Faktor of Covariance Matrix
+ 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
+ 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
+ // store fit values with tracklet
+ pTrk->SetTt(dTt); // dangerous & dirty, overwrites common fit
+ pTrk->SetTtErr(dTtErr);
+ pTrk->SetT0(dT0);
+ pTrk->SetT0Err(dT0Err);
+ pTrk->SetT0TtCov(dT0TtCov);
+ return dTt;
}
Double_t CbmTofTrackletTools::GetXdif(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit)
@@ -279,4 +286,37 @@ Double_t CbmTofTrackletTools::GetTexpected(CbmTofTracklet* pTrk, Int_t iDetId, C
return dTref;
}
+Double_t CbmTofTrackletTools::GetTexpectedError(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit, Double_t dTmean)
+{
+ Double_t dTrms = 0.;
+ Double_t dTrms2 = 0.;
+ Double_t Nref = 0;
+ Double_t dTt = 0.;
+
+ //dTt = FitTt(pTrk, iDetId);
+ dTt = pTrk->GetTt();
+
+ for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
+ if (iDetId == pTrk->GetTofDetIndex(iHit) || 0 == pTrk->GetTofDetIndex(iHit)) continue;
+ Double_t dSign = 1.;
+ if (pTrk->GetTofHitPointer(iHit)->GetZ() < pHit->GetZ()) dSign = -1;
+ Double_t dTref =
+ pTrk->GetTofHitPointer(iHit)->GetTime() - dTt * dSign * pTrk->Dist3D(pTrk->GetTofHitPointer(iHit), pHit);
+ dTref -= dTmean;
+ dTrms2 += dTref * dTref;
+ Nref++;
+ }
+ if (Nref == 0) {
+ LOG(error) << "DetId " << iDetId << ", Nref " << Nref << " sizes " << pTrk->GetNofHits();
+ return 1.E20;
+ }
+ dTrms2 /= (Double_t) Nref;
+ dTrms = TMath::Sqrt(dTrms2);
+ /*
+ LOG(info)<<"TExpErr: "<<Form("addr 0x%08x, Nhit %d, Nref %d, TM %8.2f, Rms %8.2f ",
+ iDetId, (Int_t)pTrk->GetNofHits(), (Int_t)Nref, dTmean, dTrms);
+ */
+ return dTrms;
+}
+
ClassImp(CbmTofTrackletTools)
diff --git a/reco/detectors/tof/CbmTofTrackletTools.h b/reco/detectors/tof/CbmTofTrackletTools.h
index d238d471ebcf2eaa6a52926f75a44dba4c0318f0..46c35558287d8d70205f4b93d2a1de6a559fd523 100644
--- a/reco/detectors/tof/CbmTofTrackletTools.h
+++ b/reco/detectors/tof/CbmTofTrackletTools.h
@@ -40,6 +40,7 @@ public:
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 GetTexpectedError(CbmTofTracklet* pTrk, Int_t iDetId, CbmTofHit* pHit, Double_t dTexp);
private:
ClassDef(CbmTofTrackletTools, 1);