19 #include "CbmTofDigiExp.h"
25 #include "FairEventHeader.h"
26 #include "FairFileHeader.h"
27 #include "FairGeoParSet.h"
28 #include "FairMQLogger.h"
29 #include "FairMQProgOptions.h"
30 #include "FairRootFileSink.h"
31 #include "FairRootManager.h"
32 #include "FairRunOnline.h"
33 #include "FairRuntimeDb.h"
36 #include "TClonesArray.h"
37 #include "TDirectory.h"
40 #include "TGeoManager.h"
51 #include <boost/archive/binary_iarchive.hpp>
52 #include <boost/archive/binary_oarchive.hpp>
53 #include <boost/serialization/vector.hpp>
62 using std::runtime_error::runtime_error;
93 , fTofCalDigisColl(NULL)
95 , fTofDigiMatchColl(NULL)
96 , fTofCalDigisCollOut(NULL)
97 , fTofHitsCollOut(NULL)
98 , fTofDigiMatchCollOut(NULL)
110 , fbSwapChannelSides(kFALSE)
111 , fiOutputTreeEntry(0)
137 , fhPulserTimeRawEvo()
147 , fhRpcCluDelMatPos()
150 , fhRpcCluDelMatTOff()
162 , fhRpcDTLastHits_Tot()
163 , fhRpcDTLastHits_CluSize()
165 , fhTRpcCluPosition()
176 , fhTRpcCluTOffDTLastHits()
177 , fhTRpcCluTotDTLastHits()
178 , fhTRpcCluSizeDTLastHits()
179 , fhTRpcCluMemMulDTLastHits()
219 , fdChannelDeadtime(0.)
222 , fEnableMatchPosScaling(kTRUE)
223 , fEnableAvWalk(kFALSE)
225 , fCalParFileName(
"")
226 , fOutHstFileName(
"")
227 , fOutRootFileName(
"")
230 , fRootEvent(NULL) {}
234 LOG(info) <<
"Finally close root file " <<
fOutRootFile->GetName();
253 int noChannel = fChannels.size();
254 LOG(info) <<
"Number of defined input channels: " << noChannel;
255 for (
auto const& entry : fChannels) {
256 LOG(info) <<
"Channel name: " << entry.first;
259 if (entry.first !=
"syscmd")
269 LOG(error) << e.what();
276 std::size_t pos1 = channelName.find(entry);
277 if (pos1 != std::string::npos) {
278 const vector<std::string>::const_iterator
pos =
281 LOG(info) <<
"Found " << entry <<
" in " << channelName;
282 LOG(info) <<
"Channel name " << channelName
283 <<
" found in list of allowed channel names at position "
288 LOG(info) <<
"Channel name " << channelName
289 <<
" not found in list of allowed channel names.";
290 LOG(error) <<
"Stop device.";
295 LOG(info) <<
"Init work space for CbmDeviceHitBuilderTof.";
297 fiMaxEvent = fConfig->GetValue<int64_t>(
"MaxEvent");
298 LOG(info) <<
"Max number of events to be processed: " <<
fiMaxEvent;
299 fiRunId = fConfig->GetValue<int64_t>(
"RunId");
300 fiMode = fConfig->GetValue<int64_t>(
"Mode");
301 fiPulserMode = fConfig->GetValue<int64_t>(
"PulserMode");
302 fiPulMulMin = fConfig->GetValue<uint64_t>(
"PulMulMin");
303 fiPulDetRef = fConfig->GetValue<uint64_t>(
"PulDetRef");
304 fiPulTotMin = fConfig->GetValue<uint64_t>(
"PulTotMin");
305 fiPulTotMax = fConfig->GetValue<uint64_t>(
"PulTotMax");
316 FairRunOnline* fRun =
new FairRunOnline(0);
317 rootMgr = FairRootManager::Instance();
322 if (NULL ==
fOutRootFile) LOG(fatal) <<
"could not open root file";
325 LOG(fatal) <<
"could not init Sink";
329 fDutId = fConfig->GetValue<uint64_t>(
"DutType");
330 fDutSm = fConfig->GetValue<uint64_t>(
"DutSm");
331 fDutRpc = fConfig->GetValue<uint64_t>(
"DutRpc");
333 fSelId = fConfig->GetValue<uint64_t>(
"SelType");
334 fSelSm = fConfig->GetValue<uint64_t>(
"SelSm");
335 fSelRpc = fConfig->GetValue<uint64_t>(
"SelRpc");
337 fSel2Id = fConfig->GetValue<uint64_t>(
"Sel2Type");
338 fSel2Sm = fConfig->GetValue<uint64_t>(
"Sel2Sm");
339 fSel2Rpc = fConfig->GetValue<uint64_t>(
"Sel2Rpc");
342 fiBeamRefSm = fConfig->GetValue<uint64_t>(
"BRefSm");
350 LOG(info) <<
"Init Root Output to " <<
fOutRootFile->GetName();
364 TTree* outTree =
new TTree(FairRootManager::GetTreeName(),
"/cbmout", 99);
365 LOG(info) <<
"define Tree " << outTree->GetName();
368 rootMgr->GetSink()->SetOutTree(outTree);
370 LOG(info) <<
"Initialized outTree with rootMgr at " <<
rootMgr;
384 LOG(info) <<
"Init parameter containers for CbmDeviceHitBuilderTof.";
386 FairRuntimeDb* fRtdb = FairRuntimeDb::instance();
387 if (NULL == fRtdb) LOG(error) <<
"No FairRuntimeDb found";
394 std::string parSet[NSet];
395 parSet[0] =
"CbmTofDigiPar";
396 parSet[1] =
"CbmTofDigiBdfPar";
397 parSet[2] =
"FairGeoParSet";
398 std::string Channel =
"parameters";
400 Bool_t isSimulation = kFALSE;
404 for (Int_t iSet = 0; iSet < NSet; iSet++) {
405 std::string message = parSet[iSet] +
"," + to_string(
fiRunId);
406 LOG(info) <<
"Requesting parameter container, sending message: " << message;
408 FairMQMessagePtr req(NewSimpleMessage(message));
410 FairMQMessagePtr rep(NewMessage());
412 if (Send(req, Channel) > 0) {
413 if (Receive(rep, Channel) >= 0) {
414 if (rep->GetSize() != 0) {
419 static_cast<CbmTofDigiPar*
>(tmsg.ReadObject(tmsg.GetClass()));
424 tmsg.ReadObject(tmsg.GetClass()));
426 LOG(info) <<
"Calib data file: "
440 cont =
static_cast<FairParSet*
>(tmsg.ReadObject(tmsg.GetClass()));
445 fGeoMan = (TGeoManager*) ((FairGeoParSet*) cont)
447 LOG(info) <<
"GeoMan: " <<
fGeoMan <<
" " << gGeoManager;
449 if (
k14a > iGeoVersion) {
450 LOG(error) <<
"Incompatible geometry !!!";
455 gGeoManager->Export(
"HitBuilder.geo.root");
460 LOG(warn) <<
"Parameter Set " << iSet <<
" not implemented ";
462 LOG(info) <<
"Received parameter from server:";
464 LOG(warn) <<
"Received empty reply. Parameter not available";
482 LOG(info) << Form(
"Use Dut 0x%08x, Sel 0x%08x, Sel2 0x%08x, BRef 0x%08x",
491 LOG(info) <<
"ReInit parameter containers for CbmDeviceHitBuilderTof.";
503 LOG(debug) <<
"Received message " <<
fNumMessages <<
" with " << parts.Size()
505 <<
", size0: " << parts.At(0)->GetSize();
507 std::string msgStrE(
static_cast<char*
>(parts.At(0)->GetData()),
508 (parts.At(0))->GetSize());
509 std::istringstream issE(msgStrE);
510 boost::archive::binary_iarchive inputArchiveE(issE);
520 std::string msgStr(
static_cast<char*
>(parts.At(1)->GetData()),
521 (parts.At(1))->GetSize());
522 std::istringstream iss(msgStr);
523 boost::archive::binary_iarchive inputArchive(iss);
525 std::vector<CbmTofDigiExp*> vdigi;
526 inputArchive >> vdigi;
538 LOG(debug) <<
"vdigi vector at " << vdigi.data() <<
" with size "
541 for (UInt_t iDigi = 0; iDigi < vdigi.size(); iDigi++) {
542 LOG(debug) <<
"#" << iDigi <<
" " << vdigi[iDigi]->ToString();
570 for (
int iDigi = 0; iDigi <
fiNDigiIn; iDigi++) {
572 vdigi[iDigi]->Delete();
589 LOG(info) <<
"Processed " <<
fNumMessages <<
" messages";
624 rootMgr->DeleteOldWriteoutBufferData();
629 LOG(info) <<
"File closed after " <<
fdEvent <<
" events. ";
646 const char* cmd = (
char*) (msg->GetData());
647 const char cmda[4] = {*cmd};
648 LOG(info) <<
"Handle message " << cmd <<
", " << cmd[0];
653 LOG(info) <<
"Close root file " <<
fOutRootFile->GetName();
662 if (strcmp(cmda,
"STOP")) {
678 std::this_thread::sleep_for(std::chrono::milliseconds(1000));
701 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
704 fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
707 fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
709 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
710 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
716 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
717 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] =
722 fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
723 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
724 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
725 fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
727 for (Int_t iCh = 0; iCh < iNbChan; iCh++) {
728 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
729 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
730 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
731 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
732 for (Int_t iSide = 0; iSide < nbSide; iSide++) {
733 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] =
735 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] =
737 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] =
739 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(
742 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
749 LOG(info) <<
"defaults set";
759 LOG(info) <<
"InitCalibParameter: read histos from "
767 LOG(error) <<
"InitCalibParameter: "
777 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
781 (TProfile*) gDirectory->FindObjectAny(Form(
"cl_SmT%01d_Svel", iSmType));
784 TDirectory* curdir = gDirectory;
786 gDirectory->cd(oldir->GetPath());
788 (TProfile*) hSvel->Clone();
789 gDirectory->cd(curdir->GetPath());
791 LOG(info) <<
"Svel histogram not found for module type " << iSmType;
794 for (Int_t iPar = 0; iPar < 4; iPar++) {
795 TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(
796 Form(
"cl_SmT%01d_Fpar%1d", iSmType, iPar));
797 if (NULL != hFparcur) {
798 gDirectory->cd(oldir->GetPath());
800 (TProfile*) hFparcur->Clone();
801 gDirectory->cd(curdir->GetPath());
805 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
806 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
810 if (Vscal == 0.) Vscal = 1.;
816 LOG(info) <<
"Modify " << iSmType << iSm << iRpc <<
" Svel by "
820 TH2F* htempPos_pfx = (TH2F*) gDirectory->FindObjectAny(
821 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
822 TH2F* htempTOff_pfx = (TH2F*) gDirectory->FindObjectAny(
823 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
824 TH1D* htempTot_Mean = (TH1D*) gDirectory->FindObjectAny(
825 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
826 TH1D* htempTot_Off = (TH1D*) gDirectory->FindObjectAny(
827 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
828 if (NULL != htempPos_pfx && NULL != htempTOff_pfx
829 && NULL != htempTot_Mean && NULL != htempTot_Off) {
831 Int_t iNbinTot = htempTot_Mean->GetNbinsX();
832 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
833 Double_t YMean = ((TProfile*) htempPos_pfx)
834 ->GetBinContent(iCh + 1);
836 ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
840 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0] += -dTYOff + TMean;
841 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] += +dTYOff + TMean;
843 for (Int_t iSide = 0; iSide < 2; iSide++) {
844 Double_t TotMean = htempTot_Mean->GetBinContent(
845 iCh * 2 + 1 + iSide);
846 if (0.001 < TotMean) {
847 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] *=
850 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] =
851 htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
854 if (5 == iSmType || 8 == iSmType) {
855 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] =
856 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
857 fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][1] =
859 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][1] =
860 fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][0];
863 LOG(debug) <<
"InitCalibParameter:"
864 <<
" SmT " << iSmType <<
" Sm " << iSm <<
" Rpc "
865 << iRpc <<
" Ch " << iCh
866 << Form(
": YMean %f, TMean %f", YMean, TMean) <<
" -> "
869 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][0],
870 fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][1],
873 <<
", NbinTot " << iNbinTot;
874 TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
875 Form(
"Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px",
880 TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
881 Form(
"Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px",
886 if (NULL != htempWalk0
887 && NULL != htempWalk1) {
888 LOG(debug) <<
"Initialize Walk correction for "
889 << Form(
" SmT%01d_sm%03d_rpc%03d_Ch%03d",
896 <<
"InitCalibParameter: Inconsistent Walk histograms";
898 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin] =
899 htempWalk0->GetBinContent(iBin + 1);
900 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] =
901 htempWalk1->GetBinContent(iBin + 1);
904 if (5 == iSmType || 8 == iSmType) {
905 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][1][iBin] =
906 fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][0][iBin];
912 LOG(warn) <<
" Calibration histos "
914 "cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
917 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
918 TH1D* htmpDelTof = (TH1D*) gDirectory->FindObjectAny(
919 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",
924 if (NULL == htmpDelTof) {
925 LOG(debug) <<
" Histos "
926 << Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",
934 LOG(debug) <<
" Load DelTof from histos "
935 << Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",
942 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] +=
943 htmpDelTof->GetBinContent(iBx + 1);
948 gDirectory->cd(oldir->GetPath());
949 TH1D* h1DelTof = (TH1D*) htmpDelTof->Clone(
950 Form(
"cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",
956 LOG(debug) <<
" copy histo " << h1DelTof->GetName()
957 <<
" to directory " << oldir->GetName();
959 gDirectory->cd(curdir->GetPath());
968 LOG(info) <<
"InitCalibParameter: initialization done";
981 new TH1F(
"hEvDetMul",
"Detector multiplicity of Selector; Mul", 50, 0, 50);
982 fhPulMul =
new TH1F(
"hPulMul",
"Pulser multiplicity; Mul", 110, 0, 110);
985 for (Int_t iModTyp = 0; iModTyp < 10; iModTyp++)
989 new TH2F(Form(
"hPulserTimesRaw"),
990 Form(
"Pulser Times uncorrected; Det# []; t - t0 [ns]"),
999 for (Int_t iDetSide = 0; iDetSide < iNDet * 2; iDetSide++) {
1001 Form(
"hPulserTimeRawEvo_%d", iDetSide),
1002 Form(
"Raw Pulser TimeEvolution %d; PulserEvent# ; DeltaT [ns] ",
1012 new TH2F(Form(
"hPulserTimesCor"),
1013 Form(
"Pulser Times corrected; Det# []; t - t0 [ns]"),
1022 new TH2F(Form(
"hDigiTimesRaw"),
1023 Form(
"Digi Times uncorrected; Det# []; t - t0 [ns]"),
1032 Form(
"Digi Times corrected; Det# []; t - t0 [ns]"),
1045 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
1053 Form(
"hDigiTot_SmT%01d_sm%03d_rpc%03d", iSmType, iSmId, iRpcId),
1054 Form(
"Digi Tot of Rpc #%03d in Sm %03d of type %d; digi 0; digi 1",
1066 Form(
"cl_SmT%01d_sm%03d_rpc%03d_DigiCor", iSmType, iSmId, iRpcId),
1068 "Digi Correlation of Rpc #%03d in Sm %03d of type %d; digi 0; digi 1",
1092 Double_t YSCAL = 50.;
1110 LOG(warn) <<
"No DigiPar for SmType " << Form(
"%d, 0x%08x", iS, iUCellId);
1116 new TH2F(Form(
"cl_SmT%01d_Pos", iS),
1117 Form(
"Clu position of SmType %d; Sm+Rpc# []; ypos [cm]", iS),
1125 Double_t TSumMax = 1.E3;
1128 new TH2F(Form(
"cl_SmT%01d_TOff", iS),
1129 Form(
"Clu TimeZero in SmType %d; Sm+Rpc# []; TOff [ns]", iS),
1137 TProfile* hSvelcur =
1138 (TProfile*) gDirectory->FindObjectAny(Form(
"cl_SmT%01d_Svel", iS));
1139 if (NULL == hSvelcur) {
1141 Form(
"cl_SmT%01d_Svel", iS),
1142 Form(
"Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iS),
1153 for (Int_t iPar = 0; iPar < 4; iPar++) {
1154 TProfile* hFparcur = (TProfile*) gDirectory->FindObjectAny(
1155 Form(
"cl_SmT%01d_Fpar%1d", iS, iPar));
1156 if (NULL == hFparcur) {
1157 LOG(info) <<
"Histo " << Form(
"cl_SmT%01d_Fpar%1d", iS, iPar)
1158 <<
" not found, recreate ...";
1160 Form(
"cl_SmT%01d_Fpar%1d", iS, iPar),
1161 Form(
"Clu Fpar %d in SmType %d; Sm+Rpc# []; value ", iPar, iS),
1168 fhSmCluFpar[iS][iPar] = (TProfile*) hFparcur->Clone();
1174 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
1176 new TH2F(Form(
"cl_TSmT%01d_Sel%02d_Pos", iS, iSel),
1177 Form(
"Clu position of SmType %d under Selector %02d; Sm+Rpc# "
1188 new TH2F(Form(
"cl_TSmT%01d_Sel%02d_TOff", iS, iSel),
1189 Form(
"Clu TimeZero in SmType %d under Selector %02d; Sm+Rpc# "
1200 new TH2F(Form(
"cl_TSmT%01d_Sel%02d_TRun", iS, iSel),
1201 Form(
"Clu TimeZero in SmType %d under Selector %02d; Event# "
1215 LOG(info) <<
" Define Clusterizer histos for " << iNbDet <<
" detectors ";
1239 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
1248 LOG(warn) <<
"No DigiPar for Det " << Form(
"0x%08x", iUCellId);
1251 LOG(info) <<
"DetIndx " << iDetIndx <<
", SmType " << iSmType <<
", SmId "
1252 << iSmId <<
", RpcId " << iRpcId <<
" => UniqueId "
1253 << Form(
"(0x%08x, 0x%08x)", iUniqueId, iUCellId) <<
", dx "
1265 "missing ChannelInfo for ch %d addr 0x%08x", iCh, iCCellId);
1271 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Mul", iSmType, iSmId, iRpcId),
1272 Form(
"Clu multiplicity of Rpc #%03d in Sm %03d of type %d; M []; Cnts",
1281 Form(
"cl_SmT%01d_sm%03d_rpc%03d_rate", iSmType, iSmId, iRpcId),
1282 Form(
"Clu rate of Rpc #%03d in Sm %03d of type %d; Time (s); Rate (Hz)",
1291 Form(
"cl_SmT%01d_sm%03d_rpc%03d_DTLastHits", iSmType, iSmId, iRpcId),
1292 Form(
"Clu #DeltaT to last hits of Rpc #%03d in Sm %03d of type %d; log( "
1293 "#DeltaT (ns)); counts",
1302 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Tot_DTLH", iSmType, iSmId, iRpcId),
1303 Form(
"Clu Tot of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); TOT "
1316 Form(
"cl_SmT%01d_sm%03d_rpc%03d_CluSize_DTLH", iSmType, iSmId, iRpcId),
1317 Form(
"Clu Size of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); "
1329 Double_t YSCAL = 50.;
1333 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Pos", iSmType, iSmId, iRpcId),
1335 "Clu position of Rpc #%03d in Sm %03d of type %d; Strip []; ypos [cm]",
1347 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_DelPos", iSmType, iSmId, iRpcId),
1348 Form(
"Clu position difference of Rpc #%03d in Sm %03d of type "
1349 "%d; Strip []; #Deltaypos(clu) [cm]",
1361 Form(
"cl_SmT%01d_sm%03d_rpc%03d_DelMatPos", iSmType, iSmId, iRpcId),
1362 Form(
"Matched Clu position difference of Rpc #%03d in Sm %03d of type "
1363 "%d; Strip []; #Deltaypos(mat) [cm]",
1374 Double_t TSumMax = 1.E3;
1377 Form(
"cl_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
1379 "Clu TimeZero of Rpc #%03d in Sm %03d of type %d; Strip []; TOff [ns]",
1391 Form(
"cl_SmT%01d_sm%03d_rpc%03d_DelTOff", iSmType, iSmId, iRpcId),
1392 Form(
"Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
1393 "[]; #DeltaTOff(clu) [ns]",
1405 Form(
"cl_SmT%01d_sm%03d_rpc%03d_DelMatTOff", iSmType, iSmId, iRpcId),
1406 Form(
"Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
1407 "[]; #DeltaTOff(mat) [ns]",
1419 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Trms", iSmType, iSmId, iRpcId),
1421 "Clu Time RMS of Rpc #%03d in Sm %03d of type %d; Strip []; Trms [ns]",
1433 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId),
1434 Form(
"Clu Tot of Rpc #%03d in Sm %03d of type %d; StripSide []; TOT [ns]",
1446 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Size", iSmType, iSmId, iRpcId),
1448 "Clu size of Rpc #%03d in Sm %03d of type %d; Strip []; size [strips]",
1461 Form(
"cl_SmT%01d_sm%03d_rpc%03d_AvWalk", iSmType, iSmId, iRpcId),
1462 Form(
"Walk in SmT%01d_sm%03d_rpc%03d_AvWalk", iSmType, iSmId, iRpcId),
1471 Form(
"cl_SmT%01d_sm%03d_rpc%03d_AvLnWalk", iSmType, iSmId, iRpcId),
1472 Form(
"Walk in SmT%01d_sm%03d_rpc%03d_AvLnWalk", iSmType, iSmId, iRpcId),
1483 for (Int_t iSide = 0; iSide < 2; iSide++) {
1485 new TH2D(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk",
1491 Form(
"Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk",
1514 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
1515 fhSeldT[iSel] =
new TH2F(Form(
"cl_dt_Sel%02d", iSel),
1516 Form(
"Selector time %02d; dT [ns]", iSel),
1539 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
1548 LOG(warn) <<
"No DigiPar for Det " << Form(
"0x%08x", iUCellId);
1551 LOG(debug) <<
"DetIndx " << iDetIndx <<
", SmType " << iSmType
1552 <<
", SmId " << iSmId <<
", RpcId " << iRpcId
1554 << Form(
"(0x%08x, 0x%08x)", iUniqueId, iUCellId) <<
", dx "
1574 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
1576 new TH1F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Mul",
1581 Form(
"Clu multiplicity of Rpc #%03d in Sm %03d of type %d "
1582 "under Selector %02d; M []; cnts",
1592 LOG(error) <<
" Histo not generated !";
1594 Double_t YSCAL = 50.;
1598 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Pos",
1603 Form(
"Clu position of Rpc #%03d in Sm %03d of type %d under "
1604 "Selector %02d; Strip []; ypos [cm]",
1616 Double_t TSumMax = 1.E4;
1619 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff",
1624 Form(
"Clu TimeZero of Rpc #%03d in Sm %03d of type %d under "
1625 "Selector %02d; Strip []; TOff [ns]",
1639 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot",
1644 Form(
"Clu Tot of Rpc #%03d in Sm %03d of type %d under "
1645 "Selector %02d; StripSide []; TOT [ns]",
1658 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size",
1663 Form(
"Clu size of Rpc #%03d in Sm %03d of type %d under "
1664 "Selector %02d; Strip []; size [strips]",
1678 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk",
1683 Form(
"Walk in SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk; TOT; T-TSel",
1697 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",
1702 Form(
"SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof; TRef-TSel; T-TSel",
1716 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY",
1721 Form(
"SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY; #Delta x [cm]; "
1739 for (Int_t iSide = 0; iSide < 2; iSide++) {
1741 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk",
1748 Form(
"Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk",
1765 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff_DTLH",
1770 Form(
"Clu TimeZero of Rpc #%03d in Sm %03d of type %d under "
1771 "Selector %02d; log(#DeltaT (ns)); TOff [ns]",
1784 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot_DTLH",
1789 Form(
"Clu Tot of Rpc #%03d in Sm %03d of type %d under "
1790 "Selector %02d; log(#DeltaT (ns)); TOT [ns]",
1803 new TH2F(Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size_DTLH",
1808 Form(
"Clu size of Rpc #%03d in Sm %03d of type %d under "
1809 "Selector %02d; log(#DeltaT (ns)); size [strips]",
1822 Form(
"cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_MemMul_DTLH",
1827 Form(
"Clu Memorized Multiplicity of Rpc #%03d in Sm %03d of type %d "
1828 "under Selector %02d; log(#DeltaT (ns)); TOff [ns]",
1851 TList* tHistoList(NULL);
1852 tHistoList = gROOT->GetList();
1854 TIter next(tHistoList);
1860 while ((obj = (TObject*) next())) {
1861 if (obj->InheritsFrom(TH1::Class())) {
1874 if (0 ==
fiMode)
return kTRUE;
1884 for (Int_t iDigInd = 0; iDigInd <
fiNDigiIn; iDigInd++) {
1885 CbmTofDigiExp* pDigi = &
fvDigiIn[iDigInd];
1897 Int_t iAddr = pDigi->GetAddress() &
DetMask;
1903 if (pDigi->GetTime() <
pRef->GetTime())
pRef = pDigi;
1914 fhRpcDigiTot[iDetIndx]->Fill(2 * pDigi->GetChannel() + pDigi->GetSide(),
1920 " DigiCor Histo for DetIndx %d derived from 0x%08x not found",
1922 pDigi->GetAddress());
1927 CbmTofDigiExp* pDigi2Min = NULL;
1929 for (Int_t iDigI2 = 0; iDigI2 < iNbTofDigi; iDigI2++) {
1930 CbmTofDigiExp* pDigi2 = &
fvDigiIn[iDigI2];
1932 if (0. == pDigi->GetSide() && 1. == pDigi2->GetSide()) {
1934 pDigi2->GetChannel());
1936 if (1. == pDigi->GetSide() && 0. == pDigi2->GetSide()) {
1938 pDigi->GetChannel());
1941 if (pDigi->GetSide() != pDigi2->GetSide()) {
1942 if (pDigi->GetChannel() == pDigi2->GetChannel()) {
1943 Double_t dTDif = TMath::Abs(pDigi->GetTime() - pDigi2->GetTime());
1944 if (dTDif < dTDifMin) {
1950 TMath::Abs(pDigi->GetChannel() - pDigi2->GetChannel())
1954 for (; iDigI3 < iNbTofDigi; iDigI3++) {
1955 CbmTofDigiExp* pDigi3 = &
fvDigiIn[iDigI3];
1956 if (pDigi3->GetSide() == pDigi->GetSide()
1957 && pDigi2->GetChannel() == pDigi3->GetChannel())
1960 if (iDigI3 == iNbTofDigi) {
1966 LOG(debug) << Form(
"shift channel %d%d%d%d%d and ",
1967 (Int_t) pDigi->GetType(),
1968 (Int_t) pDigi->GetSm(),
1969 (Int_t) pDigi->GetRpc(),
1970 (Int_t) pDigi->GetChannel(),
1971 (Int_t) pDigi->GetSide())
1972 << Form(
" %d%d%d%d%d ",
1973 (Int_t) pDigi2->GetType(),
1974 (Int_t) pDigi2->GetSm(),
1975 (Int_t) pDigi2->GetRpc(),
1976 (Int_t) pDigi2->GetChannel(),
1977 (Int_t) pDigi2->GetSide());
1979 if (pDigi->GetSide() == 0)
1980 pDigi2->SetAddress(pDigi->GetSm(),
1982 pDigi->GetChannel(),
1983 1 - pDigi->GetSide(),
1986 pDigi->SetAddress(pDigi2->GetSm(),
1988 pDigi2->GetChannel(),
1989 1 - pDigi2->GetSide(),
1992 LOG(debug) << Form(
"resultchannel %d%d%d%d%d and ",
1993 (Int_t) pDigi->GetType(),
1994 (Int_t) pDigi->GetSm(),
1995 (Int_t) pDigi->GetRpc(),
1996 (Int_t) pDigi->GetChannel(),
1997 (Int_t) pDigi->GetSide())
1998 << Form(
" %d%d%d%d%d ",
1999 (Int_t) pDigi2->GetType(),
2000 (Int_t) pDigi2->GetSm(),
2001 (Int_t) pDigi2->GetRpc(),
2002 (Int_t) pDigi2->GetChannel(),
2003 (Int_t) pDigi2->GetSide());
2006 CbmTofDigiExp* pDigiN =
new CbmTofDigiExp(*pDigi);
2007 pDigiN->SetAddress(pDigi->GetSm(),
2009 pDigi2->GetChannel(),
2012 pDigiN->SetTot(pDigi2->GetTot());
2015 CbmTofDigiExp* pDigiN2 =
new CbmTofDigiExp(*pDigi2);
2016 pDigiN2->SetAddress(pDigi2->GetSm(),
2018 pDigi->GetChannel(),
2021 pDigiN2->SetTot(pDigi->GetTot());
2031 if (pDigi2Min != NULL) {
2037 pDigi->GetChannel());
2041 LOG(warn) << Form(
"Invalid ChannelInfo for 0x%08x, 0x%08x",
2043 pDigi2Min->GetAddress());
2047 pDigi->GetType(), pDigi->GetSm(), pDigi->GetRpc())
2051 if (8 == pDigi->GetType() || 5 == pDigi->GetType()) {
2052 if (pDigi->GetTime() != pDigi2Min->GetTime()) {
2054 LOG(warn) <<
"Inconsistent duplicated digis in event "
2056 LOG(warn) <<
" " << pDigi->ToString();
2057 LOG(warn) <<
" " << pDigi2Min->ToString();
2059 pDigi2Min->SetTot(pDigi->GetTot());
2060 pDigi2Min->SetTime(pDigi->GetTime());
2070 for (Int_t iDigInd = 0; iDigInd <
fiNDigiIn; iDigInd++) {
2071 CbmTofDigiExp* pDigi = &
fvDigiIn[iDigInd];
2072 Int_t iAddr = pDigi->GetAddress() &
DetMask;
2074 Int_t iSide = pDigi->GetSide();
2076 pDigi->GetTime() -
pRef->GetTime());
2084 CbmTofDigiExp* pDigi;
2085 CbmTofDigiExp* pCalDigi = NULL;
2086 Int_t iDigIndCal = -1;
2088 std::map<Int_t, Double_t> mChannelDeadTime;
2090 Int_t iNbTofDigi =
fvDigiIn.size();
2091 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
2092 pDigi = (CbmTofDigiExp*) &
fvDigiIn[iDigInd];
2093 Int_t iAddr = pDigi->GetAddress();
2105 if (pDigi->GetType() == 5
2108 if (pDigi->GetSide() == 1)
2111 Bool_t bValid = kTRUE;
2112 std::map<Int_t, Double_t>::iterator it;
2113 it = mChannelDeadTime.find(iAddr);
2114 if (it != mChannelDeadTime.end()) {
2125 pCalDigi =
new ((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigiExp(*pDigi);
2127 mChannelDeadTime[iAddr] = pDigi->GetTime();
2128 if (!bValid)
continue;
2143 pCalDigi->SetTime(pCalDigi->GetTime()
2145 pCalDigi->SetTot(pCalDigi->GetTot()
2153 > pDigi->GetChannel()) {
2156 pCalDigi->GetTime() -
2159 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
2163 pCalDigi->GetTot() -
2166 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()];
2167 if (dTot < 0.001) dTot = 0.001;
2172 + pDigi->GetRpc()][pDigi->GetChannel()][pDigi->GetSide()]);
2176 Int_t iWx = (Int_t)((pCalDigi->GetTot() -
fdTOTMin) / dTotBinSize);
2177 if (0 > iWx) iWx = 0;
2180 (pCalDigi->GetTot() -
fdTOTMin) / dTotBinSize - (Double_t) iWx - 0.5;
2184 + pCalDigi->GetRpc()][pCalDigi->GetChannel()]
2185 [pCalDigi->GetSide()][iWx];
2193 + pCalDigi->GetRpc()][pCalDigi->GetChannel()]
2194 [pCalDigi->GetSide()][iWx + 1]
2198 + pCalDigi->GetRpc()][pCalDigi->GetChannel()]
2199 [pCalDigi->GetSide()][iWx]);
2208 + pCalDigi->GetRpc()][pCalDigi->GetChannel()]
2209 [pCalDigi->GetSide()][iWx - 1]
2213 + pCalDigi->GetRpc()][pCalDigi->GetChannel()]
2214 [pCalDigi->GetSide()][iWx]);
2217 pCalDigi->SetTime(pCalDigi->GetTime() - dWT);
2221 LOG(info) <<
"Skip1 Digi "
2222 <<
" Type " << pDigi->GetType() <<
" "
2225 << pDigi->GetRpc() <<
" "
2227 << pDigi->GetChannel() <<
" "
2230 if (pCalDigi->GetType() == 5
2231 || pCalDigi->GetType()
2233 CbmTofDigiExp* pCalDigi2 =
2234 new ((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigiExp(*pCalDigi);
2235 if (pCalDigi->GetSide() == 0)
2236 pCalDigi2->SetAddress(pCalDigi->GetSm(),
2238 pCalDigi->GetChannel(),
2240 pCalDigi->GetType());
2242 pCalDigi2->SetAddress(pCalDigi->GetSm(),
2244 pCalDigi->GetChannel(),
2246 pCalDigi->GetType());
2254 <<
" calibrated digis ";
2255 if (iNbTofDigi > 1) {
2259 LOG(warn) <<
"CalibRaw: Sorting not successful ";
2266 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
2268 Int_t iAddr = pDigi->GetAddress() &
DetMask;
2270 Int_t iSide = pDigi->GetSide();
2272 pDigi->GetTime() -
pRefCal->GetTime());
2282 Double_t dWeightedTime = 0.0;
2283 Double_t dWeightedPosX = 0.0;
2284 Double_t dWeightedPosY = 0.0;
2285 Double_t dWeightedPosZ = 0.0;
2286 Double_t dWeightsSum = 0.0;
2291 Int_t iNbChanInHit = 0;
2293 Int_t iLastChan = -1;
2294 Double_t dLastPosX =
2296 Double_t dLastPosY = 0.0;
2297 Double_t dLastTime = 0.0;
2299 Double_t dPosX = 0.0;
2300 Double_t dPosY = 0.0;
2301 Double_t dPosZ = 0.0;
2302 Double_t dTime = 0.0;
2303 Double_t dTimeDif = 0.0;
2304 Double_t dTotS = 0.0;
2305 Int_t fiNbSameSide = 0;
2308 gGeoManager->SetTopVolume(gGeoManager->FindVolumeFast(
"cave"));
2309 gGeoManager->CdTop();
2312 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
2315 for (Int_t iSm = 0; iSm < iNbSm; iSm++)
2316 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
2331 dWeightedTime = 0.0;
2332 dWeightedPosX = 0.0;
2333 dWeightedPosY = 0.0;
2334 dWeightedPosZ = 0.0;
2350 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
2354 if (0 ==
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc].size())
2362 1 <
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
2364 while ((
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][0])
2366 == (
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1])
2370 if (
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()
2373 <<
"SameSide Digis! on TSRC " << iSmType << iSm << iRpc
2374 << iCh <<
", Times: "
2385 << (
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][1])
2398 <<
"3 consecutive SameSide Digis! on TSRC " << iSmType
2399 << iSm << iRpc << iCh <<
", Times: "
2442 "Ev %8.0f, digis not properly time ordered, TSRCS "
2464 <<
"SameSide Erase fStor entries(d) " << iSmType
2465 <<
", SR " << iSm * iNbRpc + iRpc <<
", Ch" << iCh;
2473 if (2 >
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh]
2479 <<
"digis processing for "
2480 << Form(
" SmT %3d Sm %3d Rpc %3d Ch %3d # %3lu ",
2487 if (2 >
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh]
2490 "Leaving digi processing for TSRC %d%d%d%d, size %3lu",
2495 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size());
2506 "TSRC %d%d%d%d size %3lu ",
2511 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size())
2512 << Form(
" ChId: 0x%08x 0x%08x ", iChId, iUCellId);
2516 LOG(error) <<
"BuildHits: no geometry info! "
2517 << Form(
" %3d %3d %3d %3d 0x%08x 0x%08x ",
2535 LOG(error) << Form(
"Node at (%6.1f,%6.1f,%6.1f) : %p",
2544 CbmTofDigiExp* xDigiA =
2546 CbmTofDigiExp* xDigiB =
2549 dTimeDif = (xDigiA->GetTime() - xDigiB->GetTime());
2550 if (5 == iSmType && dTimeDif != 0.) {
2553 <<
"BuildHits: Diamond hit in " << iSm
2554 <<
" with inconsistent digits " << xDigiA->GetTime()
2555 <<
", " << xDigiB->GetTime() <<
" -> " << dTimeDif;
2561 if (1 == xDigiA->GetSide())
2571 &&
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()
2574 <<
"Hit candidate outside correlation window, check for "
2575 "better possible digis, "
2577 <<
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size();
2579 CbmTofDigiExp* xDigiC =
2581 Double_t dPosYN = 0.;
2582 Double_t dTimeDifN = 0;
2583 if (xDigiC->GetSide() == xDigiA->GetSide())
2584 dTimeDifN = xDigiC->GetTime() - xDigiB->GetTime();
2586 dTimeDifN = xDigiA->GetTime() - xDigiC->GetTime();
2588 if (1 == xDigiA->GetSide())
2595 if (TMath::Abs(dPosYN) < TMath::Abs(dPosY)) {
2597 <<
"Replace digi on side " << xDigiC->GetSide()
2598 <<
", yPosNext " << dPosYN <<
" old: " << dPosY;
2599 dTimeDif = dTimeDifN;
2601 if (xDigiC->GetSide() == xDigiA->GetSide()) {
2622 if (xDigiA->GetSide() == xDigiB->GetSide()) {
2624 <<
"Wrong combinations of digis "
2630 dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
2633 dTotS = xDigiA->GetTot() + xDigiB->GetTot();
2636 dPosX = ((Double_t)(-iNbCh / 2 + iCh) + 0.5)
2651 if (0 < iNbChanInHit) {
2652 if (iLastChan == iCh - 1) {
2665 && iLastChan == iCh - 1
2668 dWeightedTime += dTime * dTotS;
2669 dWeightedPosX += dPosX * dTotS;
2670 dWeightedPosY += dPosY * dTotS;
2671 dWeightedPosZ += dPosZ * dTotS;
2672 dWeightsSum += dTotS;
2681 <<
" Add Digi and erase fStor entries(a): NbChanInHit "
2682 << iNbChanInHit <<
", " << iSmType <<
", SR "
2683 << iSm * iNbRpc + iRpc <<
", Ch" << iCh;
2701 dWeightedTime /= dWeightsSum;
2702 dWeightedPosX /= dWeightsSum;
2703 dWeightedPosY /= dWeightsSum;
2704 dWeightedPosZ /= dWeightsSum;
2707 Double_t hitpos_local[3];
2708 hitpos_local[0] = dWeightedPosX;
2709 hitpos_local[1] = dWeightedPosY;
2710 hitpos_local[2] = dWeightedPosZ;
2714 gGeoManager->GetCurrentNode();
2715 gGeoManager->GetCurrentMatrix();
2718 gGeoManager->LocalToMaster(hitpos_local, hitpos);
2725 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
2749 if (iChm < 0) iChm = 0;
2750 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
2752 iSm, iRpc, iChm, 0, iSmType);
2771 LOG(warn) <<
"Digi refs for Hit " <<
fiNbHits <<
": "
2778 && dWeightedTime == pHitL->
GetTime()) {
2779 LOG(debug) <<
"Store Hit twice? "
2780 <<
" fiNbHits " <<
fiNbHits <<
", "
2781 << Form(
"0x%08x", iDetId);
2783 CbmTofDigiExp* pDigiC =
2786 LOG(debug) <<
" Digi " << pDigiC->ToString();
2794 CbmTofDigiExp* pDigiC =
2796 LOG(debug) <<
" DigiL " << pDigiC->ToString();
2809 Int_t(dWeightsSum * 10.));
2815 LH_store(iSmType, iSm, iRpc, iChm, pHit);
2839 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
2840 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
2850 dWeightedTime = dTime * dTotS;
2851 dWeightedPosX = dPosX * dTotS;
2852 dWeightedPosY = dPosY * dTotS;
2853 dWeightedPosZ = dPosZ * dTotS;
2854 dWeightsSum = dTotS;
2881 LOG(debug) << Form(
"1.Hit on channel %d, time: %f, PosY %f",
2887 dWeightedTime = dTime * dTotS;
2888 dWeightedPosX = dPosX * dTotS;
2889 dWeightedPosY = dPosY * dTotS;
2890 dWeightedPosZ = dPosZ * dTotS;
2891 dWeightsSum = dTotS;
2904 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2906 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2908 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2910 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin());
2928 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
2929 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].clear();
2937 <<
"=> Cluster building "
2938 <<
"from digis to hits not implemented for pads, Sm type "
2939 << iSmType <<
" Rpc " << iRpc;
2945 if (0 < iNbChanInHit) {
2957 dWeightedTime /= dWeightsSum;
2958 dWeightedPosX /= dWeightsSum;
2959 dWeightedPosY /= dWeightsSum;
2960 dWeightedPosZ /= dWeightsSum;
2963 Double_t hitpos_local[3] = {3 * 0.};
2964 hitpos_local[0] = dWeightedPosX;
2965 hitpos_local[1] = dWeightedPosY;
2966 hitpos_local[2] = dWeightedPosZ;
2970 gGeoManager->GetCurrentNode();
2971 gGeoManager->GetCurrentMatrix();
2975 gGeoManager->LocalToMaster(hitpos_local, hitpos);
2983 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
2996 if (iChm < 0) iChm = 0;
2997 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
3014 LOG(warn) <<
"Digi refs for Hit " <<
fiNbHits <<
": "
3020 && dWeightedTime == pHitL->
GetTime())
3021 LOG(debug) <<
"Store Hit twice? "
3022 <<
" fiNbHits " <<
fiNbHits <<
", "
3023 << Form(
"0x%08x", iDetId);
3034 Int_t(dWeightsSum * 10.));
3041 LH_store(iSmType, iSm, iRpc, iChm, pHit);
3063 fvdX[iSmType][iRpc].push_back(dWeightedPosX);
3064 fvdY[iSmType][iRpc].push_back(dWeightedPosY);
3090 if (
fvLastHits[iSmType][iSm][iRpc][iChm].size() == 0)
3091 fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
3093 Double_t dLastTime = pHit->
GetTime();
3094 if (dLastTime >=
fvLastHits[iSmType][iSm][iRpc][iChm].back()->GetTime()) {
3095 fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit);
3097 " Store LH from Ev %8.0f for TSRC %d%d%d%d, size %lu, addr 0x%08x, "
3107 dLastTime -
fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime());
3113 std::list<CbmTofHit*>::iterator it;
3114 for (it =
fvLastHits[iSmType][iSm][iRpc][iChm].begin();
3115 it !=
fvLastHits[iSmType][iSm][iRpc][iChm].end();
3117 if ((*it)->GetTime() > dLastTime)
break;
3118 fvLastHits[iSmType][iSm][iRpc][iChm].insert(--it, pHit);
3119 Double_t deltaTime = dLastTime - (*it)->GetTime();
3120 LOG(debug) << Form(
"Hit inserted into LH from Ev %8.0f for TSRC "
3121 "%d%d%d%d, size %lu, addr 0x%08x, delta time %f ",
3131 Double_t deltaTime =
3132 dLastTime -
fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime();
3133 LOG(debug) << Form(
"first LH from Ev %8.0f for TSRC %d%d%d%d, size "
3134 "%lu, addr 0x%08x, delta time %f ",
3143 if (deltaTime == 0.) {
3147 fvLastHits[iSmType][iSm][iRpc][iChm].push_front(pHit);
3156 LOG(error) << Form(
"Inconsistent LH Smtype size %lu, %d ",
3162 LOG(error) << Form(
"Inconsistent LH Sm size %lu, %d T %d",
3169 LOG(error) << Form(
"Inconsistent LH Rpc size %lu, %d TS %d%d ",
3175 if ((Int_t)
fvLastHits[iSmType][iSm][iRpc].size()
3178 "Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
3185 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
3191 "Inconsistent address for Ev %8.0f in list of size %lu for "
3192 "TSRC %d%d%d%d: 0x%08x, time %f",
3200 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
3205 LOG(debug) << Form(
"LH check passed for event %8.0f",
fdEvent);
3210 LOG(error) << Form(
"Inconsistent LH Smtype size %lu, %d ",
3215 LOG(error) << Form(
"Inconsistent LH Sm size %lu, %d T %d",
3222 LOG(error) << Form(
"Inconsistent LH Rpc size %lu, %d TS %d%d ",
3228 if ((Int_t)
fvLastHits[iSmType][iSm][iRpc].size()
3231 "Inconsistent LH RpcChannel size %lu, %d TSR %d%d%d ",
3238 while (
fvLastHits[iSmType][iSm][iRpc][iCh].size() > 0) {
3244 "Inconsistent address for Ev %8.0f in list of size %lu for "
3245 "TSRC %d%d%d%d: 0x%08x, time %f",
3253 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
3254 fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete();
3255 fvLastHits[iSmType][iSm][iRpc][iCh].pop_front();
3260 LOG(info) << Form(
"LH cleaning done after %8.0f events",
fdEvent);
3270 Double_t dLastTotS) {
3276 Int_t iCh = iLastChan + 1;
3277 LOG(debug) << Form(
"Inspect channel TSRC %d%d%d%d at time %f, pos %f, size ",
3284 <<
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size();
3285 if (iCh == iNbCh)
return kFALSE;
3286 if (0 ==
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size())
3290 if (1 <
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
3291 Bool_t AddedHit = kFALSE;
3294 < (Int_t)
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size() - 1;
3296 if (AddedHit)
break;
3301 < (Int_t)
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].size()) {
3304 if ((
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][i1])->GetSide()
3305 == (
fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh][i2])
3311 CbmTofDigiExp* xDigiA =
3313 CbmTofDigiExp* xDigiB =
3315 Double_t dTime = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
3320 gGeoManager->FindNode(
3323 Double_t dTimeDif = xDigiA->GetTime() - xDigiB->GetTime();
3324 Double_t dPosY = 0.;
3325 if (1 == xDigiA->GetSide())
3331 if (TMath::Abs(dPosY - dLastPosY)
3334 Double_t dNClHits = (Double_t)(
vDigiIndRef.size() / 2);
3337 Double_t dTotS = xDigiA->GetTot() + xDigiB->GetTot();
3338 Double_t dNewTotS = (dLastTotS + dTotS);
3339 dLastPosX = (dLastPosX * dLastTotS + dPosX * dTotS) / dNewTotS;
3340 dLastPosY = (dLastPosY * dLastTotS + dPosY * dTotS) / dNewTotS;
3341 dLastTime = (dLastTime * dLastTotS + dTime * dTotS) / dNewTotS;
3342 dLastTotS = dNewTotS;
3344 Int_t Ind1 =
fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i1];
3345 Int_t Ind2 =
fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh][i2];
3350 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + i1);
3352 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + i1);
3354 std::vector<int>::iterator it;
3355 it = find(
fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin(),
3358 if (it !=
fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].end()) {
3360 it -
fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin();
3361 LOG(debug) <<
"Found i2 " << i2 <<
" with Ind2 " << Ind2
3362 <<
" at position " << ipos;
3364 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + ipos);
3366 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc][iCh].begin() + ipos);
3368 LOG(error) <<
" Did not find i2 " << i2 <<
" with Ind2 " << Ind2;
3372 if (iCh != (iNbCh - 1)
3381 LOG(debug) <<
"Added Strip " << iCh <<
" to cluster of size "
3392 Double_t hitpos_local[3] = {3 * 0.};
3393 hitpos_local[0] = dLastPosX;
3394 hitpos_local[1] = dLastPosY;
3395 hitpos_local[2] = 0.;
3398 TGeoNode* cNode = gGeoManager->GetCurrentNode();
3399 if (NULL == cNode) {
3400 LOG(error) << Form(
"Node at (%6.1f,%6.1f,%6.1f) : %p",
3409 gGeoManager->GetCurrentMatrix();
3410 gGeoManager->LocalToMaster(hitpos_local, hitpos);
3411 TVector3 hitPos(hitpos[0], hitpos[1], hitpos[2]);
3412 TVector3 hitPosErr(0.5, 0.5, 0.5);
3414 if (iChm < 0) iChm = 0;
3415 if (iChm > iNbCh - 1) iChm = iNbCh - 1;
3437 Int_t(dLastTotS * 10.));
3441 LH_store(iSmType, iSm, iRpc, iChm, pHit);
3460 double wx = 1. - par[4] * TMath::Power(xx + par[5], 2);
3461 double xboxe = par[0] * 0.25
3462 * (1. + TMath::Erf((xx + par[1] - par[3]) / par[2]))
3463 * (1. + TMath::Erf((-xx + par[1] + par[3]) / par[2]));
3469 h1 = (TH1*) gROOT->FindObjectAny(hname);
3470 if (NULL != h1) {
fit_ybox(h1, 0.); }
3474 Double_t* fpar = NULL;
3480 Double_t* fpar = NULL) {
3481 TAxis* xaxis = h1->GetXaxis();
3482 Double_t Ymin = xaxis->GetXmin();
3483 Double_t Ymax = xaxis->GetXmax();
3484 TF1* f1 =
new TF1(
"YBox",
f1_xboxe, Ymin, Ymax, 6);
3485 Double_t yini = (h1->GetMaximum() + h1->GetMinimum()) * 0.5;
3486 if (ysize == 0.) ysize = Ymax * 0.8;
3487 f1->SetParameters(yini, ysize * 0.5, 1., 0., 0., 0.);
3489 f1->SetParLimits(2, 0.2, 3.);
3490 f1->SetParLimits(3, -4., 4.);
3493 for (Int_t
i = 0;
i < 4;
i++)
3495 for (Int_t
i = 0;
i < 4;
i++)
3496 f1->SetParameter(2 +
i, fp[
i]);
3497 LOG(debug) <<
"Ini Fpar for " << h1->GetName() <<
" with "
3498 << Form(
" %6.3f %6.3f %6.3f %6.3f ", fp[0], fp[1], fp[2], fp[3]);
3501 h1->Fit(
"YBox",
"Q");
3505 res[9] = f1->GetChisquare();
3507 for (
int i = 0;
i < 6;
i++) {
3508 res[
i] = f1->GetParameter(
i);
3509 err[
i] = f1->GetParError(
i);
3512 LOG(debug) <<
"YBox Fit of " << h1->GetName()
3513 <<
" ended with chi2 = " << res[9]
3514 << Form(
", strip length %7.2f +/- %5.2f, position resolution "
3515 "%7.2f +/- %5.2f at y_cen = %7.2f +/- %5.2f",
3527 <<
" geometrically known modules ";
3532 LOG(info) <<
"Digi Parameter container contains " << iNrOfCells <<
" cells.";
3534 for (Int_t icell = 0; icell < iNrOfCells; ++icell) {
3550 LOG(debug) <<
"-I- InitPar " << icell <<
" Id: " << Form(
"0x%08x", cellId)
3551 <<
" " << cell <<
" tmcs: " << smtype <<
" " << smodule <<
" "
3552 << module <<
" " << cell <<
" x=" << Form(
"%6.2f",
x)
3553 <<
" y=" << Form(
"%6.2f",
y) <<
" z=" << Form(
"%6.2f", z)
3554 <<
" dx=" << dx <<
" dy=" << dy;
3557 gGeoManager->FindNode(
3560 if (NULL == fNode) {
3561 LOG(error) << Form(
"Node at (%6.1f,%6.1f,%6.1f) : %p",
3570 TGeoHMatrix* cMatrix = gGeoManager->GetCurrentMatrix();
3581 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
3586 LOG(info) <<
" DetIndx " << iDetIndx <<
"(" << iNbDet <<
"), SmType "
3587 << iSmType <<
", SmId " << iSmId <<
", RpcId " << iRpcId
3588 <<
" => UniqueId " << Form(
"0x%08x ", iUniqueId)
3589 << Form(
" Svel %6.6f, DeadStrips 0x%08x ",
3602 for (Int_t
i = 0;
i < 2;
i++)
3616 fvdX.resize(iNbSmTypes);
3617 fvdY.resize(iNbSmTypes);
3624 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
3631 fvdX[iSmType].resize(iNbRpc);
3632 fvdY[iSmType].resize(iNbRpc);
3633 fvdDifX[iSmType].resize(iNbRpc);
3634 fvdDifY[iSmType].resize(iNbRpc);
3636 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
3639 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
3644 LOG(warn) <<
"LoadGeometry: StoreDigi without channels "
3645 << Form(
"SmTy %3d, Sm %3d, NbRpc %3d, Rpc, %3d ",
3651 fStorDigiExp[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
3652 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
3653 fvLastHits[iSmType][iSm][iRpc].resize(iNbChan);
3663 fvdX.resize(iNbSmTypes);
3664 fvdY.resize(iNbSmTypes);
3668 for (Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
3671 fStorDigi[iSmType].resize(iNbSm * iNbRpc);
3675 fvdX[iSmType].resize(iNbRpc);
3676 fvdY[iSmType].resize(iNbRpc);
3677 fvdDifX[iSmType].resize(iNbRpc);
3678 fvdDifY[iSmType].resize(iNbRpc);
3680 for (Int_t iSm = 0; iSm < iNbSm; iSm++) {
3681 for (Int_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
3683 fStorDigi[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
3684 fStorDigiInd[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
3697 CbmTofDigiExp* pDigi;
3700 for (Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
3718 > pDigi->GetChannel()) {
3721 + pDigi->GetRpc()][pDigi->GetChannel()]
3725 + pDigi->GetRpc()][pDigi->GetChannel()]
3726 .push_back(iDigInd);
3728 LOG(info) <<
"Skip2 Digi "
3729 <<
" Type " << pDigi->GetType() <<
" "
3732 << pDigi->GetRpc() <<
" "
3734 << pDigi->GetChannel() <<
" "
3743 for (Int_t iHit = 0; iHit <
fiNbHits; iHit++) {
3746 "Found Hit %d, addr 0x%08x, X %6.2f, Y %6.2f Z %6.2f T %6.2f CLS %d",
3756 std::vector<CbmTofHit*> vhit;
3758 for (Int_t iHit = 0; iHit <
fiNbHits; iHit++) {
3764 std::stringstream ossE;
3765 boost::archive::binary_oarchive oaE(ossE);
3767 std::string* strMsgE =
new std::string(ossE.str());
3769 std::stringstream oss;
3770 boost::archive::binary_oarchive oa(oss);
3772 std::string* strMsg =
new std::string(oss.str());
3775 parts.AddPart(NewMessage(
3776 const_cast<char*
>(strMsgE->c_str()),
3778 [](
void*,
void*
object) { delete static_cast<std::string*>(object); },
3781 parts.AddPart(NewMessage(
3782 const_cast<char*
>(strMsg->c_str()),
3784 [](
void*,
void*
object) { delete static_cast<std::string*>(object); },
3787 if (Send(parts,
"tofhits")) {
3788 LOG(error) <<
"Problem sending data ";
3801 gGeoManager->CdTop();
3803 if (0 < iNbTofHits) {
3805 Double_t dTTrig[
iNSel];
3807 Double_t ddXdZ[
iNSel];
3808 Double_t ddYdZ[
iNSel];
3809 Double_t dSel2dXdYMin[
iNSel];
3811 Int_t iBeamRefMul = 0;
3812 Int_t iBeamAddRefMul = 0;
3816 LOG(debug) <<
"FillHistos() for " <<
iNSel <<
" triggers"
3830 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
3842 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
3844 if (NULL == pHit)
continue;
3851 std::map<UInt_t, UInt_t>::iterator it =
fDetIdIndexMap.find(iDetId);
3854 Int_t iDetIndx = it->second;
3866 &&
fvLastHits[iSmType][iSm][iRpc][iCh].size() == 0)
3867 LOG(error) << Form(
" <E> hit not stored in memory for TSRC %d%d%d%d",
3874 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
3883 LOG(warn) << Form(
"Invalid time ordering in ev %8.0f in list of "
3884 "size %lu for TSRC %d%d%d%d: Delta t %f ",
3893 while (
fvLastHits[iSmType][iSm][iRpc][iCh].size() > 2.
3894 ||
fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()
3898 <<
" pop from list size "
3899 <<
fvLastHits[iSmType][iSm][iRpc][iCh].size()
3900 << Form(
" outdated hits for ev %8.0f in TSRC %d%d%d%d",
3907 " with tHit - tLast %f ",
3909 -
fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime())
3915 "Inconsistent address in list of size %lu for TSRC %d%d%d%d: "
3923 fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime());
3924 fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete();
3925 fvLastHits[iSmType][iSm][iRpc][iCh].pop_front();
3930 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
3933 Double_t dTotSum = 0.;
3934 for (Int_t iLink = 0; iLink < digiMatch->
GetNofLinks();
3938 Int_t iDigInd1 = (digiMatch->
GetLink(iLink + 1)).GetIndex();
3939 CbmTofDigiExp* pDig0 =
3941 CbmTofDigiExp* pDig1 =
3943 dTotSum += pDig0->GetTot() + pDig1->GetTot();
3946 std::list<CbmTofHit*>::iterator itL =
3950 iH <
fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1;
3954 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()));
3956 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()),
3959 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()), dTotSum);
3967 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
3969 if (NULL == pHit)
continue;
3978 Double_t TotSum = 0.;
3979 for (Int_t iLink = 0; iLink < digiMatch->
GetNofLinks();
3983 if (iDigInd0 < fTofCalDigisColl->GetEntries()) {
3984 CbmTofDigiExp* pDig0 =
3986 TotSum += pDig0->GetTot();
4001 LOG(debug) <<
"FillHistos: BRefMul: " << iBeamRefMul <<
", "
4003 if (iBeamRefMul == 0)
4008 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
4009 BSel[iSel] = kFALSE;
4017 dSel2dXdYMin[iSel] = 1.E300;
4026 for (Int_t iRpc = iR0; iRpc < iRl; iRpc++)
4028 LOG(debug) <<
"Selector 0: DutMul "
4030 << iDutMul <<
", BRefMul " << iBeamRefMul
4031 <<
" TRef: " <<
dTRef <<
", BeamAddRefMul "
4036 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
4038 if (NULL == pHit)
continue;
4045 if (pHit->
GetTime() < dTTrig[iSel]) {
4046 dTTrig[iSel] = pHit->
GetTime();
4053 "Found selector 0 with mul %d from 0x%08x at %f ",
4066 for (Int_t iRpc = iR0; iRpc < iRl; iRpc++)
4068 LOG(debug) <<
"FillHistos(): selector 1: RefMul "
4070 << iRefMul <<
", BRefMul " << iBeamRefMul;
4074 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
4076 if (NULL == pHit)
continue;
4080 if (pHit->
GetTime() < dTTrig[iSel]) {
4081 dTTrig[iSel] = pHit->
GetTime();
4088 "Found selector 1 with mul %d from 0x%08x at %f ",
4096 LOG(info) <<
"FillHistos: selection not implemented " << iSel;
4103 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
4105 BSel[iSel] = kFALSE;
4108 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
4110 if (NULL == pHit)
continue;
4113 Double_t dzscal = 1.;
4115 dzscal = pHit->
GetZ() / pTrig[iSel]->
GetZ();
4116 Double_t dSEl2dXdz = (pHit->
GetX() - pTrig[iSel]->
GetX())
4117 / (pHit->
GetZ() - pTrig[iSel]->
GetZ());
4118 Double_t dSEl2dYdz = (pHit->
GetY() - pTrig[iSel]->
GetY())
4119 / (pHit->
GetZ() - pTrig[iSel]->
GetZ());
4123 pHit->
GetX() - dzscal * pTrig[iSel]->
GetX(), 2.)
4125 pHit->
GetY() - dzscal * pTrig[iSel]->
GetY(), 2.))
4128 Double_t dX2Y2 = TMath::Sqrt(dSEl2dXdz * dSEl2dXdz
4129 + dSEl2dYdz * dSEl2dYdz);
4130 if (dX2Y2 < dSel2dXdYMin[iSel]) {
4131 ddXdZ[iSel] = dSEl2dXdz;
4132 ddYdZ[iSel] = dSEl2dYdz;
4133 dSel2dXdYMin[iSel] = dX2Y2;
4142 UInt_t uTriggerPattern = 1;
4144 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4146 uTriggerPattern |= (0x1 << (iSel * 3
4151 for (Int_t iSel = 0; iSel <
iNSel; iSel++) {
4154 for (UInt_t uChannel = 0; uChannel < 16; uChannel++) {
4155 if (uTriggerPattern & (0x1 << uChannel)) {
4164 for (Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) {
4166 if (NULL == pHit)
continue;
4169 std::map<UInt_t, UInt_t>::iterator it =
fDetIdIndexMap.find(iDetId);
4172 Int_t iDetIndx = it->second;
4179 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4193 if (iBeamRefMul == 0)
break;
4199 LOG(error) <<
"Invalid Channel Pointer for ChId "
4200 << Form(
" 0x%08x ", iChId) <<
", Ch " << iCh;
4204 gGeoManager->FindNode(
4213 hitpos[0] = pHit->
GetX();
4214 hitpos[1] = pHit->
GetY();
4215 hitpos[2] = pHit->
GetZ();
4216 Double_t hitpos_local[3];
4218 gGeoManager->GetCurrentNode();
4219 gGeoManager->MasterToLocal(hitpos, hitpos_local);
4228 (Double_t) iCh, hitpos_local[1]);
4232 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4249 LOG(error) <<
" Inconsistent DigiMatches for Hitind " << iHitInd
4261 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4265 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
4267 std::list<CbmTofHit*>::iterator itL =
4271 iH <
fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1;
4275 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()),
4281 Double_t TotSum = 0.;
4282 for (Int_t iLink = 0; iLink < digiMatch->
GetNofLinks();
4286 if (iDigInd0 < fTofCalDigisColl->GetEntries()) {
4287 CbmTofDigiExp* pDig0 =
4289 TotSum += pDig0->GetTot();
4292 Double_t dMeanTimeSquared = 0.;
4293 Double_t dNstrips = 0.;
4295 Double_t dDelTof = 0.;
4296 Double_t dTcor[
iNSel];
4297 Double_t dTTcor[
iNSel];
4298 Double_t dZsign[
iNSel];
4299 Double_t dzscal = 1.;
4302 for (Int_t iLink = 0; iLink < digiMatch->
GetNofLinks();
4307 (digiMatch->
GetLink(iLink + 1)).GetIndex();
4310 if (iDigInd0 < fTofCalDigisColl->GetEntries()
4312 CbmTofDigiExp* pDig0 =
4314 CbmTofDigiExp* pDig1 =
4316 if ((Int_t) pDig0->GetType() != iSmType) {
4317 LOG(error) << Form(
" Wrong Digi SmType for Tofhit %d in iDetIndx "
4318 "%d, Ch %d with %3.0f strips at Indx %d, %d",
4334 pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
4336 pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
4338 Int_t iCh0 = pDig0->GetChannel();
4339 Int_t iCh1 = pDig1->GetChannel();
4340 Int_t iS0 = pDig0->GetSide();
4341 Int_t iS1 = pDig1->GetSide();
4342 if (iCh0 != iCh1 || iS0 == iS1) {
4344 " MT2 for Tofhit %d in iDetIndx %d, Ch %d from %3.0f strips: ",
4349 << Form(
" Dig0: Ind %d, Ch %d, Side %d, T: %6.1f ",
4354 << Form(
" Dig1: Ind %d, Ch %d, Side %d, T: %6.1f ",
4364 " Wrong Digi for Tofhit %d in iDetIndx %d, Ch %d at Indx %d, %d "
4365 "from %3.0f strips: %d, %d, %d, %d",
4384 dMeanTimeSquared += TMath::Power(
4385 0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->
GetTime(), 2);
4389 TMath::Log10(0.5 * (pDig0->GetTot() + pDig1->GetTot())),
4390 0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->
GetTime());
4392 Double_t dTotWeigth = (pDig0->GetTot() + pDig1->GetTot()) / TotSum;
4393 Double_t dCorWeigth = 1. - dTotWeigth;
4396 pDig0->GetChannel(),
4398 * (0.5 * (pDig0->GetTime() + pDig1->GetTime())
4401 Double_t dDelPos = 0.5 * (pDig0->GetTime() - pDig1->GetTime())
4403 if (0 == pDig0->GetSide()) dDelPos *= -1.;
4405 pDig0->GetChannel(), dCorWeigth * (dDelPos - hitpos_local[1]));
4411 - (1. - 2. * pDig0->GetSide()) * hitpos_local[1]
4417 - (1. - 2. * pDig1->GetSide()) * hitpos_local[1]
4424 - (1. - 2. * pDig0->GetSide()) * hitpos_local[1]
4430 - (1. - 2. * pDig1->GetSide()) * hitpos_local[1]
4440 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4442 if (NULL == pHit || NULL == pTrig[iSel]) {
4443 LOG(info) <<
" invalid pHit, iSel " << iSel <<
", iDetIndx "
4450 pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
4452 pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
4453 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
4455 std::list<CbmTofHit*>::iterator itL =
4459 iH <
fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1;
4463 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()),
4466 TMath::Log10(pHit->
GetTime() - (*itL)->GetTime()),
4472 dzscal = pHit->
GetZ() / pTrig[iSel]->
GetZ();
4474 pHit->
GetX() - dzscal * pTrig[iSel]->
GetX(),
4475 pHit->
GetY() - dzscal * pTrig[iSel]->
GetY());
4477 if (pHit->
GetZ() < pTrig[iSel]->
GetZ()) dZsign[iSel] = -1.;
4483 TMath::Power(pHit->
GetX() - dzscal * pTrig[iSel]->
GetX(),
4486 pHit->
GetY() - dzscal * pTrig[iSel]->
GetY(), 2.))
4492 - (pTrig[iSel]->
GetX()
4494 * (pHit->
GetZ() - (pTrig[iSel]->
GetZ()))),
4498 - (pTrig[iSel]->
GetY()
4500 * (pHit->
GetZ() - (pTrig[iSel]->
GetZ()))),
4507 && TMath::Abs(
dTRef - dTTrig[iSel])
4525 if (iBx < 0) iBx = 0;
4529 dTentry - ((Double_t) iBx) * dBinWidth;
4531 dDTentry < 0 ? iBx1 = iBx - 1 : iBx1 = iBx + 1;
4532 Double_t w0 = 1. - TMath::Abs(dDTentry) / dBinWidth;
4533 Double_t w1 = 1. - w0;
4534 if (iBx1 < 0) iBx1 = 0;
4537 fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx][iSel] * w0
4538 +
fvCPDelTof[iSmType][iSm * iNbRpc + iRpc][iBx1][iSel]
4547 dTTcor[iSel] = dDelTof * dZsign[iSel];
4548 dTcor[iSel] = pHit->
GetTime() - dDelTof - dTTrig[iSel];
4553 " TRpcCluWalk for Ev %d, Link %d(%d), Sel %d, TSR %d%d%d, "
4554 "Ch %d,%d, S %d,%d T %f, DelTof %6.1f, W-ent: %6.0f,%6.0f",
4574 << Form(
" Inconsistent walk histograms -> debugging "
4575 "necessary ... for %d, %d, %d, %d, %d, %d, %d ",
4618 (Double_t)(iSm * iNbRpc + iRpc), dTcor[iSel]);
4624 hitpos[0] = pHit->
GetX() - dzscal * pTrig[iSel]->
GetX()
4626 hitpos[1] = pHit->
GetY() - dzscal * pTrig[iSel]->
GetY()
4628 hitpos[2] = pHit->
GetZ();
4629 gGeoManager->MasterToLocal(
4630 hitpos, hitpos_local);
4635 (pHit->
GetTime() - dTTrig[iSel]) - dTTcor[iSel]);
4642 LOG(error) <<
"FillHistos: invalid digi index " << iDetIndx
4643 <<
" digi0,1" << iDigInd0 <<
", " << iDigInd1
4651 Double_t dVar = dMeanTimeSquared / (dNstrips - 1);
4653 Double_t dTrms = TMath::Sqrt(dVar);
4654 LOG(debug) << Form(
" Trms for Tofhit %d in iDetIndx %d, Ch %d from "
4655 "%3.0f strips: %6.3f ns",
4671 if (TMath::Abs(hitpos_local[1])
4675 fhSmCluTOff[iSmType]->Fill((Double_t)(iSm * iNbRpc + iRpc),
4678 for (Int_t iSel = 0; iSel <
iNSel; iSel++)
4687 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
4689 std::list<CbmTofHit*>::iterator itL =
4693 iH <
fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1;
4705 if (
fvLastHits[iSmType][iSm][iRpc][iCh].size()
4707 std::list<CbmTofHit*>::iterator itL =
4710 for (Int_t iH = 0; iH < 1; iH++) {
4713 Double_t dTsinceLast = pHit->
GetTime() - (*itL)->GetTime();
4715 LOG(error) << Form(
"Invalid Time since last hit on channel "
4716 "TSRC %d%d%d%d: %f > %f",
4725 TMath::Log10(dTsinceLast), pHit->
GetTime() - dTTrig[iSel]);
4727 TMath::Log10(dTsinceLast),
4728 fvLastHits[iSmType][iSm][iRpc][iCh].size() - 1);
4743 const Double_t Tlim = 0.5;
4749 LOG(debug) <<
"Incomplete or distorted pulser event " <<
iNPulserFound
4750 <<
" with " <<
fiNDigiIn <<
" digis instead of "
4762 for (Int_t iDigi = 0; iDigi < (Int_t)
fiNDigiIn; iDigi++) {
4763 Int_t iCh =
fvDigiIn[iDigi].GetChannel();
4764 if (iCh != 0 && iCh != 31)
continue;
4768 if (iDet == iDet0) {
4769 Int_t iSide =
fvDigiIn[iDigi].GetSide();
4779 "Pulser reference %d not found in pulser event %d of mul %d, return",
4787 LOG(debug) <<
fiNDigiIn <<
" pulser digis with ref in "
4788 << Form(
"0x%08x at %d with deltaT %12.3f",
4794 for (Int_t iDigi = 0; iDigi <
fiNDigiIn; iDigi++) {
4795 Int_t iCh =
fvDigiIn[iDigi].GetChannel();
4804 Int_t iSide =
fvDigiIn[iDigi].GetSide();
4809 if (
fvDigiIn[iDigi].GetType() != 5) {
4810 if (iCh != 0 && iCh != 31)
continue;
4815 LOG(debug) <<
"Found PulHit of Type 5 in Ch " << iCh;
4818 if (iCh > 39) iSide = 1;
4819 if (iCh != 0 && iCh != 40)
continue;
4824 if (
fvDigiIn[iDigi].GetType() == 8)
4825 if (
fvDigiIn[iDigi].GetRpc() != 7)
continue;
4826 if (iCh != 0 && iCh != 31)
continue;
4832 if (iDet == iDet0 && 1 == iSide) {
4838 if (TMath::Abs(
fvDigiIn[iDigi].GetTime()
4845 if (TMath::Abs(
fvDigiIn[iDigi].GetTime()
4848 LOG(info) << Form(
"Unexpected Pulser Time for event %d, pulser %d: "
4849 "%15.1f instead %15.1f, TDif: %10.1f != %10.1f",
4866 (Double_t)(
fvDigiIn[iDigi].GetTime()));
4874 LOG(info) <<
"Unexpected pulser offset at ev " <<
fdEvent
4876 << Form(
", addr 0x%08x", iAddr) <<
", side " << iSide
4896 for (Int_t iDet = 0; iDet < (Int_t)
fvPulserTimes.size(); iDet++) {
4897 for (Int_t iSide = 0; iSide < 2; iSide++) {
4898 if (iDet == iDet0 && iSide == 1)
continue;
4900 Double_t Tmean = 0.;
4901 std::list<Double_t>::iterator it;
4909 LOG(debug) <<
"New pulser offset at ev " <<
fdEvent <<
", pulcnt "
4911 << iSide <<
": " << Tmean <<
" ( "
4919 for (Int_t iDigi = 0; iDigi <
fiNDigiIn; iDigi++) {
4921 Int_t iCh =
fvDigiIn[iDigi].GetChannel();
4924 Int_t iSide =
fvDigiIn[iDigi].GetSide();
4928 if (
fvDigiIn[iDigi].GetType() != 5) {
4929 if (iCh != 0 && iCh != 31)
continue;
4932 if (iCh > 39) iSide = 1;
4934 if (iCh != 0 && iCh != 40)
continue;
4938 if (
fvDigiIn[iDigi].GetType() == 8)
4939 if (
fvDigiIn[iDigi].GetRpc() != 7)
continue;
4940 if (iCh != 0 && iCh != 31)
continue;
4955 for (Int_t iDigi = 0; iDigi <
fiNDigiIn; iDigi++) {
4958 Int_t iSide =
fvDigiIn[iDigi].GetSide();
4961 if (5 ==
fvDigiIn[iDigi].GetType()) {
4963 Int_t iCh =
fvDigiIn[iDigi].GetChannel();
4964 if (iCh > 4) iSide = 1;
4968 const Int_t iRefAddr = 0x00078006;