CbmRoot
Loading...
Searching...
No Matches
CbmTofCalibrator.cxx
Go to the documentation of this file.
1/* Copyright (C) 2020-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Norbert Herrmann, Florian Uhlig [committer] */
4
11// CBMroot classes and includes
12#include "CbmTofCalibrator.h"
13
14#include "CbmEvent.h"
15#include "CbmMatch.h"
16#include "CbmTofAddress.h" // in cbmdata/tof
17#include "CbmTofCell.h" // in tof/TofData
18#include "CbmTofClusterizersDef.h"
19#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
20#include "CbmTofDigiBdfPar.h" // in tof/TofParam
21#include "CbmTofDigiPar.h" // in tof/TofParam
22#include "CbmTofGeoHandler.h" // in tof/TofTools
23#include "CbmTofHit.h"
24#include "CbmTofTracklet.h"
25
26// FAIR classes and includes
27#include "FairRunAna.h"
28#include "FairRuntimeDb.h"
29
30#include <Logger.h>
31
32// ROOT Classes and includes
33#include "TF1.h"
34
35#include <TClonesArray.h>
36#include <TDirectory.h>
37#include <TFile.h>
38#include <TFitResult.h>
39#include <TGeoManager.h>
40#include <TH1.h>
41#include <TH2.h>
42#include <TH3.h>
43#include <TMinuit.h>
44#include <TProfile.h>
45#include <TROOT.h>
46
47const int MaxShift = 10; // number of 50ps bins
48const double cLight = 29.9792; // in cm/ns
49const TString cBadChannelFile = "TofBadChannels.txt";
50const double dValidEdge = 1.0; // FIXME: constants in code
51
52std::vector<TH1*> hSvel;
53
54CbmTofCalibrator::CbmTofCalibrator()
55 : fDigiMan(NULL)
56 , fTofClusterizer(NULL)
57 , fTofFindTracks(NULL)
58 , fTrackletTools(NULL)
59 , fDigiPar(NULL)
60 , fDigiBdfPar(NULL)
61 , fTofDigiMatchColl(NULL)
62 , fhCalR0(NULL)
63 , fhCalDX0(NULL)
64 , fhCalDY0(NULL)
65 , fhCalCounterDt(NULL)
66 , fhCalCounterDy(NULL)
67 , fhCalChannelDt(NULL)
68 , fhCalChannelDy(NULL)
69 , fhCalTot()
70 , fhCalPosition()
71 , fhCalPos()
72 , fhCalTOff()
73 , fhCalTofOff()
74 , fhCalDelPos()
75 , fhCalDelTOff()
76 , fhCalCluTrms()
77 , fhCalCluSize()
78 , fhCalWalk()
79 , fhCalDtWalk()
80 , fhCorPos()
81 , fhCorTOff()
82 , fhCorTot()
83 , fhCorTotOff()
84 , fhCorSvel()
85 , fhCorWalk()
86 , fDetIdIndexMap()
87 , fhDoubletDt()
88 , fhDoubletDd()
89 , fhDoubletV()
90{
91}
92
93CbmTofCalibrator::~CbmTofCalibrator() {}
94
95InitStatus CbmTofCalibrator::Init()
96{
97
98 FairRootManager* fManager = FairRootManager::Instance();
99 // Get access to TofCalDigis
100 fTofCalDigiVec = fManager->InitObjectAs<std::vector<CbmTofDigi> const*>("TofCalDigi");
101 //dynamic_cast<std::vector<CbmTofDigi> const*>(fManager->GetObject("TofCalDigi"));
102 if (NULL == fTofCalDigiVec) LOG(warning) << "No access to TofCalDigis!";
103
104 // check for availability of digis
105 fDigiMan = CbmDigiManager::Instance();
106 if (NULL == fDigiMan) {
107 LOG(warning) << "No CbmDigiManager";
108 return kFATAL;
109 }
110 fDigiMan->Init();
111
112 if (!fDigiMan->IsPresent(ECbmModuleId::kTof)) {
113 LOG(warn) << "CbmTofCalibrator: No digi input!";
114 }
115
116 // check for access to current calibration file
117 fTofClusterizer = CbmTofEventClusterizer::Instance();
118 if (NULL == fTofClusterizer) {
119 LOG(warn) << "CbmTofCalibrator: no access to active calibration";
120 }
121 else {
122 TString CalParFileName = fTofClusterizer->GetCalParFileName();
123 LOG(info) << "Current calibration file: " << CalParFileName;
124 }
125 fTofFindTracks = CbmTofFindTracks::Instance();
126 if (NULL == fTofFindTracks) {
127 LOG(warn) << "CbmTofCalibrator: no access to tof tracker ";
128 }
129
130 fEvtHeader = (FairEventHeader*) fManager->GetObject("EventHeader.");
131 if (NULL == fEvtHeader) {
132 LOG(warning) << "CbmTofCalibrator::RegisterInputs => Could not get EvtHeader Object";
133 return kFATAL;
134 }
135 // Get Access to MatchCollection
136 fTofDigiMatchColl = (TClonesArray*) fManager->GetObject("TofHitCalDigiMatch");
137 if (NULL == fTofDigiMatchColl) {
138 LOG(warn) << "No branch TofHitCalDigiMatch found, try TofHitDigiMatch";
139 fTofDigiMatchColl = (TClonesArray*) fManager->GetObject("TofHitDigiMatch");
140 }
141
142 if (NULL == fTofDigiMatchColl) {
143 LOG(error) << "CbmTofCalibrator: no access to DigiMatch ";
144 return kFATAL;
145 }
146 LOG(info) << Form("Got DigiMatchColl %s at %p", fTofDigiMatchColl->GetName(), fTofDigiMatchColl);
147
148 // Get Parameters
149 if (!InitParameters()) {
150 LOG(error) << "CbmTofCalibrator: No parameters!";
151 }
152 // generate deviation histograms
153 if (!CreateCalHist()) {
154 LOG(error) << "CbmTofCalibrator: No histograms!";
155 return kFATAL;
156 }
157
158 fTofId = new CbmTofDetectorId_v21a();
159
160 fTrackletTools = new CbmTofTrackletTools(); // initialize tools
161
162 TString shcmd = Form("rm %s ", cBadChannelFile.Data());
163 gSystem->Exec(shcmd.Data());
164
165 shcmd = Form("touch %s ", cBadChannelFile.Data());
166 gSystem->Exec(shcmd.Data());
167
168 LOG(info) << "TofCalibrator initialized successfully";
169 return kSUCCESS;
170}
171
172Bool_t CbmTofCalibrator::InitParameters()
173{
174 LOG(info) << "InitParameters: get par pointers from RTDB";
175 // Get Base Container
176 FairRun* ana = FairRun::Instance();
177 FairRuntimeDb* rtdb = ana->GetRuntimeDb();
178
179 fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
180 if (NULL == fDigiPar) return kFALSE;
181
182 fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
183 if (NULL == fDigiBdfPar) return kFALSE;
184
185 return kTRUE;
186}
187
188Bool_t CbmTofCalibrator::CreateCalHist()
189{
190 TDirectory* oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
191 gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
192
193 // detector related distributions
194 Int_t iNbDet = fDigiBdfPar->GetNbDet();
195 LOG(info) << "Define Calibrator histos for " << iNbDet << " detectors "
196 << "in directory " << gDirectory->GetName() << ", old " << oldir->GetName();
197
198 fhCalR0 = new TH1D("hCalR0", "Tracklet distance to nominal vertex; R_0 [cm]", 100, 0., 0.5);
199 fhCalDX0 = new TH1D("hCalDX0", "Tracklet distance to nominal vertex; #DeltaX_0 [cm]", 100, -1.5, 1.5);
200 fhCalDY0 = new TH1D("hCalDY0", "Tracklet distance to nominal vertex; #DeltaY_0 [cm]", 100, -1.5, 1.5);
201
202 fhCalCounterDt = new TH1D("hCalCounterDt", "CalibCounterDt ; #Deltat [ns]", 100, -0.2, 0.2);
203 fhCalCounterDy = new TH1D("hCalCounterDy", "CalibCounterDy ; #Deltay [cm]", 100, -0.8, 0.8);
204 fhCalChannelDt = new TH1D("hCalChannelDt", "CalibChannelDt ; #Deltat [ns]", 100, -0.25, 0.25);
205 fhCalChannelDy = new TH1D("hCalChannelDy", "CalibChannelDy ; #Deltat [ns]", 100, -1.8, 1.8);
206
207 fhCalTot.resize(iNbDet);
208 fhCalPosition.resize(iNbDet);
209 fhCalPos.resize(iNbDet);
210 fhCalTOff.resize(iNbDet);
211 fhCalTofOff.resize(iNbDet);
212 fhCalDelPos.resize(iNbDet);
213 fhCalDelTOff.resize(iNbDet);
214 fhCalCluTrms.resize(iNbDet);
215 fhCalCluSize.resize(iNbDet);
216 fhCalWalk.resize(iNbDet);
217 fhCalDtWalk.resize(iNbDet);
218 fhCalWalkAv.resize(iNbDet);
219
220 fhCalXYTOff.resize(iNbDet);
221 fhCalXYTot.resize(iNbDet);
222 fhCalTotYWalk.resize(iNbDet);
223 fhCalTotYTOff.resize(iNbDet);
224
225 fDetIdIndexMap.clear();
226
227 for (Int_t iDetIndx = 0; iDetIndx < iNbDet; iDetIndx++) {
228 Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
229 fDetIdIndexMap[iUniqueId] = iDetIndx;
230
231 Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
232 Int_t iSmId = CbmTofAddress::GetSmId(iUniqueId);
233 Int_t iRpcId = CbmTofAddress::GetRpcId(iUniqueId);
234 Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSmId, iRpcId, 0, 0, iSmType);
235 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iUCellId);
236 if (NULL == fChannelInfo) {
237 LOG(warning) << "No DigiPar for Det " << Form("0x%08x", iUCellId);
238 continue;
239 }
240 Double_t YMAX = TMath::Max(2., fChannelInfo->GetSizey()) * 0.75;
241
242 Double_t TotMax = 20.;
243 fhCalTot[iDetIndx] = new TH2F(
244 Form("cal_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId),
245 Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; Strip []; Tot [a.u.]", iRpcId, iSmId, iSmType),
246 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, 0., TotMax);
247
248 fhCalXYTot[iDetIndx] =
249 new TH3F(Form("calXY_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId),
250 Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; Strip []; y [cm]; Tot [a.u.]", iRpcId, iSmId, iSmType),
251 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, 2 * fDigiBdfPar->GetNbChan(iSmType, iRpcId), 20, -YMAX,
252 YMAX, 100, 0., TotMax);
253
254 fhCalPosition[iDetIndx] =
255 new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_Position", iSmType, iSmId, iRpcId),
256 Form("Clu position of Rpc #%03d in Sm %03d of type %d; Strip []; ypos [cm]", iRpcId, iSmId, iSmType),
257 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YMAX, YMAX);
258
259 Double_t YDMAX = 5;
260 fhCalPos[iDetIndx] =
261 new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_Pos", iSmType, iSmId, iRpcId),
262 Form("Clu pos deviation of Rpc #%03d in Sm %03d of type %d; "
263 "Strip []; ypos [cm]",
264 iRpcId, iSmId, iSmType),
265 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -YDMAX, YDMAX);
266
267 Double_t TSumMax = 2.;
268 //if(iSmType == 5) TSumMax *= 2.; // enlarge for diamond / beamcounter
269 fhCalTOff[iDetIndx] = new TH2F(
270 Form("cal_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
271 Form("Clu T0 deviation of Rpc #%03d in Sm %03d of type %d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType),
272 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 200, -TSumMax, TSumMax);
273
274 fhCalXYTOff[iDetIndx] = new TH3F(
275 Form("calXY_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId),
276 Form("XY T0 deviation of Rpc #%03d in Sm %03d of type %d; Strip []; y[cm]; TOff [ns]", iRpcId, iSmId, iSmType),
277 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 100, -YMAX, YMAX, 100,
278 -TSumMax * 0.5, TSumMax * 0.5);
279
280 TSumMax = 20.;
281 fhCalTofOff[iDetIndx] = new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_TofOff", iSmType, iSmId, iRpcId),
282 Form("Clu T0 deviation of Rpc #%03d in Sm %03d of type %d; "
283 "Strip []; TofOff [ns]",
284 iRpcId, iSmId, iSmType),
285 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0,
286 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 199, -TSumMax, TSumMax);
287
288 fhCalDelPos[iDetIndx] =
289 new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_DelPos", iSmType, iSmId, iRpcId),
290 Form("Clu position difference of Rpc #%03d in Sm %03d of type "
291 "%d; Strip []; #Deltaypos(clu) [cm]",
292 iRpcId, iSmId, iSmType),
293 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -10., 10.);
294
295 fhCalDelTOff[iDetIndx] =
296 new TH2F(Form("cal_SmT%01d_sm%03d_rpc%03d_DelTOff", iSmType, iSmId, iRpcId),
297 Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip "
298 "[]; #DeltaTOff(clu) [ns]",
299 iRpcId, iSmId, iSmType),
300 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, -0.6, 0.6);
301
302 fhCalCluTrms[iDetIndx] =
303 new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Trms", iSmType, iSmId, iRpcId),
304 Form("Clu Time RMS of Rpc #%03d in Sm %03d of type %d; Strip []; Trms [ns]", iRpcId, iSmId, iSmType),
305 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 99, 0., 0.5);
306
307 fhCalCluSize[iDetIndx] =
308 new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Size", iSmType, iSmId, iRpcId),
309 Form("Clu size of Rpc #%03d in Sm %03d of type %d; Strip []; size [strips]", iRpcId, iSmId, iSmType),
310 fDigiBdfPar->GetNbChan(iSmType, iRpcId), 0, fDigiBdfPar->GetNbChan(iSmType, iRpcId), 16, 0.5, 16.5);
311
312 TSumMax = 0.8;
313 fhCalWalkAv[iDetIndx] =
314 new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_WalkAv", iSmType, iSmId, iRpcId),
315 Form("Walk in SmT%01d_sm%03d_rpc%03d_WalkAv; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId),
316 nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
317
318 fhCalWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
319 fhCalDtWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
320 fhCalTotYWalk[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
321 fhCalTotYTOff[iDetIndx].resize(fDigiBdfPar->GetNbChan(iSmType, iRpcId));
322
323 for (Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpcId); iCh++) {
324 fhCalWalk[iDetIndx][iCh].resize(2);
325 fhCalDtWalk[iDetIndx][iCh].resize(2);
326 fhCalTotYWalk[iDetIndx][iCh].resize(2);
327 fhCalTotYTOff[iDetIndx][iCh].resize(2);
328 for (Int_t iSide = 0; iSide < 2; iSide++) {
329 fhCalWalk[iDetIndx][iCh][iSide] =
330 new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
331 Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId,
332 iRpcId, iCh, iSide),
333 nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
334 fhCalDtWalk[iDetIndx][iCh][iSide] =
335 new TH2D(Form("cal_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_DtWalk", iSmType, iSmId, iRpcId, iCh, iSide),
336 Form("SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_DtWalk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId,
337 iCh, iSide),
338 nbClWalkBinX, 0., TotMax, nbClWalkBinY, -TSumMax, TSumMax);
339
340 fhCalTotYWalk[iDetIndx][iCh][iSide] =
341 new TH3D(Form("calTotY_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide),
342 Form("YWalk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; y [cm]; #DeltaT [ns]", iSmType,
343 iSmId, iRpcId, iCh, iSide),
344 nbClWalkBinX, 0., TotMax, 20, -YMAX, YMAX, nbClWalkBinY, -0.4, 0.4);
345
346 fhCalTotYTOff[iDetIndx][iCh][iSide] =
347 new TH3D(Form("calTotY_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_TOff", iSmType, iSmId, iRpcId, iCh, iSide),
348 Form("YTot in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_TOff; Tot [a.u.]; y [cm]; #DeltaT [ns]", iSmType,
349 iSmId, iRpcId, iCh, iSide),
350 nbClWalkBinX, 0., TotMax, 20, -YMAX, YMAX, nbClWalkBinY, -TSumMax * 0.5, TSumMax * 0.5);
351 }
352 }
353 }
354
355 return kTRUE;
356}
357
358static Int_t NevtH = 0;
359void CbmTofCalibrator::FillCalHist(CbmTofHit* pRef, Int_t iOpt, CbmEvent* tEvent)
360{
361 if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! for Event " << tEvent;
362
363 NevtH++;
364 if (NevtH == 1) LOG(info) << "FillCalHist: look for beam counter at " << pRef;
365
366 double dTAv = 0.;
367 if (pRef != NULL)
368 dTAv = pRef->GetTime();
369 else {
370 if (-1 < fTofClusterizer->GetBeamAddr()) { // look for beam counter
371 if (NevtH == 1) LOG(info) << Form("FillCalHist: look for beam counter 0x%08x", fTofClusterizer->GetBeamAddr());
372 Int_t iHit = 0;
373 for (; iHit < (int) fTofClusterizer->GetNbHits(); iHit++) {
374 CbmTofHit* pHit = fTofClusterizer->GetHitPointer(iHit);
375 if (NULL == pHit) continue;
376 if (pHit->GetAddress() == fTofClusterizer->GetBeamAddr()) {
377 pRef = pHit;
378 dTAv = pRef->GetTime();
379 if (NevtH == 1) LOG(info) << "FillCalHist: found beam counter";
380 break;
381 }
382 }
383 if (iHit == fTofClusterizer->GetNbHits()) return; // beam counter not present
384 }
385 else {
386 double dNAv = 0.;
387 for (Int_t iHit = 0; iHit < (int) fTofClusterizer->GetNbHits(); iHit++) {
388 LOG(debug) << "GetHitPointer for " << iHit;
389 CbmTofHit* pHit = fTofClusterizer->GetHitPointer(iHit);
390 if (NULL == pHit) continue;
391 dTAv += pHit->GetTime();
392 dNAv++;
393 }
394 if (dNAv == 0) return;
395 dTAv /= dNAv;
396 }
397 }
398
399 for (Int_t iHit = 0; iHit < (int) fTofClusterizer->GetNbHits(); iHit++) {
400 CbmTofHit* pHit = fTofClusterizer->GetHitPointer(iHit);
401 if (NULL == pHit) continue;
402
403 Int_t iDetId = (pHit->GetAddress() & DetMask);
404 Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
405 Int_t iSm = CbmTofAddress::GetSmId(iDetId);
406 Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
407
408 std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
409 if (it == fDetIdIndexMap.end()) {
410 LOG(info) << "detector TSR " << iSmType << iSm << iRpc << " not found in detector map";
411 continue; // continue for invalid detector index
412 }
413 Int_t iDetIndx = it->second; //fDetIdIndexMap[iDetId];
414
415 Int_t iChId = pHit->GetAddress();
416 Int_t iCh = CbmTofAddress::GetChannelId(iChId);
417
418 // fill CalDigi Tots
419 CbmMatch* digiMatch = fTofClusterizer->GetMatchPointer(iHit);
420 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // loop over digis
421 CbmLink L0 = digiMatch->GetLink(iLink); //vDigish.at(ivDigInd);
422 UInt_t iDigInd0 = L0.GetIndex();
423 UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex(); //vDigish.at(ivDigInd+1);
424 //LOG(debug1)<<" " << iDigInd0<<", "<<iDigInd1;
425 if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
426 const CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
427 const CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
428 if ((Int_t) pDig0->GetType() != iSmType) {
429 LOG(error) << Form(" Wrong Digi SmType %d - %d for Tofhit %d in iDetIndx "
430 "%d, Ch %d with %d strips at Indx %d, %d",
431 iSmType, (Int_t) pDig0->GetType(), iHit, iDetIndx, iCh,
432 fDigiBdfPar->GetNbChan(iSmType, iRpc), iDigInd0, iDigInd1);
433
434 for (Int_t iL = 0; iL < digiMatch->GetNofLinks(); iL++) { // loop over digis
435 int idx = (digiMatch->GetLink(iL)).GetIndex();
436 const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
437 LOG(info) << Form(" Digi %d, idx %d, TSRCS %d%d%d%02d%d ", iL, idx, (int) pDx->GetType(),
438 (int) pDx->GetSm(), (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
439 }
440
441 for (Int_t idx = 0; idx < (int) fTofCalDigiVec->size(); idx++) { // loop over digis
442 const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
443 LOG(info) << Form(" AllDigi %d, TSRCS %d%d%d%02d%d ", idx, (int) pDx->GetType(), (int) pDx->GetSm(),
444 (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
445 }
446
447 LOG(fatal) << " fhRpcCluTot: Digi 0 " << iDigInd0 << ": Ch " << pDig0->GetChannel() << ", Side "
448 << pDig0->GetSide() << ", StripSide " << (Double_t) iCh * 2. + pDig0->GetSide() << " Digi 1 "
449 << iDigInd1 << ": Ch " << pDig1->GetChannel() << ", Side " << pDig1->GetSide() << ", StripSide "
450 << (Double_t) iCh * 2. + pDig1->GetSide() << ", Tot0 " << pDig0->GetTot() << ", Tot1 "
451 << pDig1->GetTot();
452 } // end of type mismatch condition
453 fhCalTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
454 fhCalTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
455 }
456 else {
457 LOG(warn) << "Invalid CalDigi index " << iDigInd0 << ", " << iDigInd1 << " for size "
458 << fTofCalDigiVec->size();
459 }
460 }
461
462 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
463 if (NULL == fChannelInfo) {
464 LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
465 continue;
466 }
467 //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
468 Double_t hitpos[3];
469 hitpos[0] = pHit->GetX();
470 hitpos[1] = pHit->GetY();
471 hitpos[2] = pHit->GetZ();
472 Double_t hlocal_p[3];
473 //TGeoNode* cNode=
474 //gGeoManager->GetCurrentNode();
475 //gGeoManager->MasterToLocal(hitpos, hlocal_p);
476 fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
477
478 /*
479 LOG(info)<<Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f ",
480 NevtH,iSmType,iSm,iRpc,iCh,hlocal_p[1],pHit->GetTime() - dTAv, pRef->GetTime());
481 */
482 fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
483 /*
484 if(pHit != pRef )
485 fhCalTofOff[iDetIndx]->Fill((Double_t) iCh,
486 pHit->GetTime() - pRef->GetTime()); // ignore any dependence
487 */
488 if (iOpt < 100) {
489 fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
490 fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
491 fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
492 }
493 else {
494 double dTexp = pHit->GetR() * fTofFindTracks->GetTtTarg();
495 fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
496 fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
497 fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
498 }
499 }
500
501 return;
502}
503
504void CbmTofCalibrator::FillHitCalHist(CbmTofHit* pRef, const Int_t iOpt, CbmEvent* tEvent, TClonesArray* tTofHitsColl)
505{
506 if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! for Event " << tEvent;
507 if (NULL == tTofHitsColl) LOG(fatal) << "No access to TofHits! for Event " << tEvent;
508
509 NevtH++;
510 if (NevtH == 1)
511 LOG(warn) << "HitCalHist at " << NevtH << ", Opt = " << iOpt << ", pRef " << pRef << ", "
512 << Form("0x%08x", fTofClusterizer->GetBeamAddr());
513
514 double dTAv = 0.;
515 if (pRef != NULL) {
516 dTAv = pRef->GetTime();
517 //LOG(debug) << "HitCalHist got Tbeam " << dTAv << ", Opt = " << iOpt;
518 }
519 else {
520 if (-1 < fTofClusterizer->GetBeamAddr()) { // look for beam counter
521 if (NevtH == 1) LOG(warn) << Form("FillHitCalHist: look for beam counter 0x%08x", fTofClusterizer->GetBeamAddr());
522 size_t iHit = 0;
523 //for (; iHit < (int) fTofFindTracks->GetNbHits(); iHit++) {
524 for (; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
525 //CbmTofHit* pHit = fTofFindTracks->GetHitPointer(iHit);
526 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
527 CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
528 if (NULL == pHit) continue;
529 if ((pHit->GetAddress() & DetMask) == (fTofClusterizer->GetBeamAddr() & DetMask)) {
530 if (pRef != NULL) {
531 if (pHit->GetTime() < dTAv) pRef = pHit;
532 }
533 else {
534 pRef = pHit;
535 }
536 dTAv = pRef->GetTime();
537 if (NevtH == 1) LOG(warn) << "FillHitCalHist: found beam counter";
538 break;
539 }
540 }
541 if (iHit == tEvent->GetNofData(ECbmDataType::kTofHit)) return; // beam counter not present
542 }
543 else {
544 double dNAv = 0.;
545 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
546 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
547 CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
548 if (NULL == pHit) continue;
549 dTAv += pHit->GetTime();
550 dNAv++;
551 }
552 if (dNAv == 0) return;
553 dTAv /= dNAv;
554 }
555 }
556
557 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
558 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
559 CbmTofHit* pHit = (CbmTofHit*) tTofHitsColl->At(iHitIndex);
560 if (NULL == pHit) continue;
561
562 Int_t iDetId = (pHit->GetAddress() & DetMask);
563 Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
564 Int_t iSm = CbmTofAddress::GetSmId(iDetId);
565 Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
566
567 std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
568 if (it == fDetIdIndexMap.end()) {
569 LOG(info) << "detector TSR " << iSmType << iSm << iRpc << " not found in detector map";
570 continue; // continue for invalid detector index
571 }
572 Int_t iDetIndx = it->second; //fDetIdIndexMap[iDetId];
573
574 Int_t iChId = pHit->GetAddress();
575 Int_t iCh = CbmTofAddress::GetChannelId(iChId);
576 double dTexp = 0.;
577
578 /*
579 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
580 if (NULL == fChannelInfo) {
581 LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
582 continue;
583 }
584 */
585
586 Double_t hitpos[3];
587 hitpos[0] = pHit->GetX();
588 hitpos[1] = pHit->GetY();
589 hitpos[2] = pHit->GetZ();
590 Double_t hlocal_p[3];
591
592 fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
593
594 LOG(debug) << Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f ", NevtH, iSmType, iSm, iRpc, iCh,
595 hlocal_p[1], pHit->GetTime() - dTAv, pRef->GetTime());
596
597 fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
598
599 if (iOpt < 100) {
600 LOG(debug) << "Fill " << fhCalTofOff[iDetIndx]->GetName();
601 if (pHit != pRef) {
602 fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
603 fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv);
604 fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv);
605 }
606 }
607 else {
608 if (iSmType != 5) {
609 if (NULL != pRef) {
610 if ((pHit->GetAddress() & DetMask) != (pRef->GetAddress() & DetMask)) {
611 dTexp = pHit->Dist3D(pRef) / cLight; // Edge expected for speed of light,
612 if (pHit->GetZ() < pRef->GetZ()) dTexp *= -1;
613
614 LOG(debug) << Form(" Ev %d: TSRC %d%d%d%2d, y %6.2f, tof %6.2f, Tref %12.1f, Texp %12.1f ", NevtH, iSmType,
615 iSm, iRpc, iCh, hlocal_p[1], pHit->GetTime() - dTAv, pRef->GetTime(), dTexp);
616 }
617 }
618 else {
619 dTexp = pHit->GetR() / cLight; // Edge expected for speed of light, fTofFindTracks->GetTtTarg();
620 }
621 /*
622 LOG(debug)<<"Fill " << fhCalTofOff[iDetIndx]->GetName()<<" " << iHit <<": " << iCh << ", "
623 << pHit->GetTime() - dTAv <<" - " << dTexp;
624 */
625 if (dTexp != 0.) {
626 fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
627 fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1], pHit->GetTime() - dTAv - dTexp);
628 fhCalTofOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTAv - dTexp);
629 }
630 }
631 }
632
633 // fill CalDigi Tots and internal walks
634 if (dTexp != 0.) {
635 CbmMatch* digiMatch = fTofClusterizer->GetMatchIndexPointer(iHitIndex);
636 if (NULL == digiMatch) LOG(fatal) << "No digiMatch for Hit " << iHit << ", TSR " << iSmType << iSm << iRpc;
637 if (digiMatch->GetNofLinks() < fTofClusterizer->GetCluSizeMin()) continue;
638 double dMeanTimeSquared = 0.;
639 double dTotSum = 0.;
640 double dYSum = 0.;
641 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // first loop over digis
642 CbmLink L0 = digiMatch->GetLink(iLink);
643 UInt_t iDigInd0 = L0.GetIndex();
644 UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
645 //LOG(debug1)<<" " << iDigInd0<<", "<<iDigInd1;
646 if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
647 const CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
648 const CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
649 if ((Int_t) pDig0->GetType() != iSmType) {
650 LOG(error) << Form(" Wrong Digi SmType %d - %d for Tofhit %lu in iDetIndx "
651 "%d, Ch %d with %d strips at Indx %d, %d",
652 iSmType, (Int_t) pDig0->GetType(), iHit, iDetIndx, iCh,
653 fDigiBdfPar->GetNbChan(iSmType, iRpc), iDigInd0, iDigInd1);
654
655 for (Int_t iL = 0; iL < digiMatch->GetNofLinks(); iL++) { // loop over digis
656 int idx = (digiMatch->GetLink(iL)).GetIndex();
657 const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
658 LOG(info) << Form(" Digi %d, idx %d, TSRCS %d%d%d%02d%d ", iL, idx, (int) pDx->GetType(),
659 (int) pDx->GetSm(), (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
660 }
661
662 for (Int_t idx = 0; idx < (int) fTofCalDigiVec->size(); idx++) { // loop over digis
663 const CbmTofDigi* pDx = &(fTofCalDigiVec->at(idx));
664 LOG(info) << Form(" AllDigi %d, TSRCS %d%d%d%02d%d ", idx, (int) pDx->GetType(), (int) pDx->GetSm(),
665 (int) pDx->GetRpc(), (int) pDx->GetChannel(), (int) pDx->GetSide());
666 }
667
668 LOG(fatal) << " fhRpcCluTot: Digi 0 " << iDigInd0 << ": Ch " << pDig0->GetChannel() << ", Side "
669 << pDig0->GetSide() << ", StripSide " << (Double_t) iCh * 2. + pDig0->GetSide() << " Digi 1 "
670 << iDigInd1 << ": Ch " << pDig1->GetChannel() << ", Side " << pDig1->GetSide() << ", StripSide "
671 << (Double_t) iCh * 2. + pDig1->GetSide() << ", Tot0 " << pDig0->GetTot() << ", Tot1 "
672 << pDig1->GetTot();
673 } // end of type mismatch condition
674 fhCalTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), pDig0->GetTot());
675 fhCalTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), pDig1->GetTot());
676 if (digiMatch->GetNofLinks() > 2) {
677 double dTotStrip = pDig0->GetTot() + pDig1->GetTot();
678 dTotSum += dTotStrip;
679 double dYdigi = 0.5 * (pDig0->GetTime() - pDig1->GetTime()) * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
680 if (pDig0->GetSide() == 1) {
681 dYSum += (pDig0->GetTime() - pDig1->GetTime()) * dTotStrip; // weighted mean
682 }
683 else {
684 dYSum += (pDig1->GetTime() - pDig0->GetTime()) * dTotStrip;
685 dYdigi *= -1;
686 }
687 fhCalXYTot[iDetIndx]->Fill(pDig0->GetChannel() * 2. + pDig0->GetSide(), dYdigi, pDig0->GetTot());
688 fhCalXYTot[iDetIndx]->Fill(pDig1->GetChannel() * 2. + pDig1->GetSide(), dYdigi, pDig1->GetTot());
689 }
690 }
691 } // first loop over digis
692 if (digiMatch->GetNofLinks() > 2) {
693 double dYAv = dYSum / dTotSum;
694 double dYLoc = 0.;
695 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // second loop over digis
696 CbmLink L0 = digiMatch->GetLink(iLink);
697 UInt_t iDigInd0 = L0.GetIndex();
698 UInt_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
699 if (iDigInd0 < fTofCalDigiVec->size() && iDigInd1 < fTofCalDigiVec->size()) {
700 const CbmTofDigi* pDig0 = &(fTofCalDigiVec->at(iDigInd0));
701 const CbmTofDigi* pDig1 = &(fTofCalDigiVec->at(iDigInd1));
702 double dTotStrip = pDig0->GetTot() + pDig1->GetTot();
703 double dTotWeight = 1 - dTotStrip / dTotSum;
704 //dTotWeight *= 0.5;
705 dTotWeight = 1.;
706 double dDeltaT = 0.5 * (pDig0->GetTime() + pDig1->GetTime()) - pHit->GetTime();
707 //if(NULL != fTofFindTracks) dDeltaT += fTofFindTracks->GetTOff(iDetId); //obsolete ??
708 fhCalDelTOff[iDetIndx]->Fill(pDig0->GetChannel(), dDeltaT * dTotWeight);
709 // TBD: position deviation in counter reference frame
710 Double_t dDelPos = 0.5 * (pDig0->GetTime() - pDig1->GetTime()) * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
711 if (0 == pDig0->GetSide()) dDelPos *= -1.;
712 fhCalDelPos[iDetIndx]->Fill(pDig0->GetChannel(), (dDelPos - fTofClusterizer->GetLocalY(pHit)));
713
714 int iCh0 = pDig0->GetChannel();
715 int iSide0 = pDig0->GetSide();
716 int iCh1 = pDig1->GetChannel();
717 int iSide1 = pDig1->GetSide();
718 if (iCh0 != iCh1) LOG(fatal) << "Inconsistent TofHit - FixIt";
719 fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(), dDeltaT * dTotWeight);
720 fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(), dDeltaT * dTotWeight);
721 dMeanTimeSquared += dDeltaT * dDeltaT;
722 if (pDig0->GetSide() == 1) {
723 dYLoc = pDig0->GetTime() - pDig1->GetTime();
724 }
725 else {
726 dYLoc = pDig1->GetTime() - pDig0->GetTime();
727 }
728 dYLoc *= 0.5;
729 fhCalTotYWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(),
730 dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), dDeltaT);
731 fhCalTotYWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(),
732 dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), dDeltaT);
733
734 fhCalTotYTOff[iDetIndx][iCh0][iSide0]->Fill(
735 pDig0->GetTot(), dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), pHit->GetTime() - dTAv - dTexp);
736 fhCalTotYTOff[iDetIndx][iCh1][iSide1]->Fill(
737 pDig1->GetTot(), dYLoc * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc), pHit->GetTime() - dTAv - dTexp);
738
739 Double_t dDeltaDT = dYLoc - dYAv;
740 fhCalPos[iDetIndx]->Fill(iCh0, dDeltaDT * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
741
742 fhCalDtWalk[iDetIndx][iCh0][iSide0]->Fill(pDig0->GetTot(), (2. * pDig0->GetSide() - 1) * dDeltaDT);
743 fhCalDtWalk[iDetIndx][iCh1][iSide1]->Fill(pDig1->GetTot(), (2. * pDig1->GetSide() - 1) * dDeltaDT);
744 }
745 } // second loop over digis
746 } // if (digiMatch->GetNofLinks()>2)
747
748 fhCalCluSize[iDetIndx]->Fill((Double_t) iCh, digiMatch->GetNofLinks() / 2);
749 if (digiMatch->GetNofLinks() > 2) {
750 double dVar = dMeanTimeSquared / (digiMatch->GetNofLinks() / 2 - 1);
751 double dTrms = TMath::Sqrt(dVar);
752 fhCalCluTrms[iDetIndx]->Fill((Double_t) iCh, dTrms);
753 }
754 } // iOpt<100
755 }
756 return;
757} // -- FillHitCalHist
758
759
760static Int_t NevtT = 0;
761void CbmTofCalibrator::FillCalHist(CbmTofTracklet* pTrk, Int_t iOpt, CbmEvent* tEvent)
762{
763 NevtT++;
764 if (NULL == fTofCalDigiVec) LOG(fatal) << "No access to TofCalDigis! Opt " << iOpt << ", Event" << tEvent;
765 // fill deviation histograms on walk level
766 /*
767 LOG(info)<<"Calibrator " << NevtT <<": Tt " << pTrk->GetTt()
768 <<", R0: " << pTrk->GetR0() <<", bB "
769 << pTrk->ContainsAddr(CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 5));
770 */
771 if (pTrk->GetTt() < 0) return; // take tracks with positive velocity only
772 if (pTrk->GetNofHits() < fTofFindTracks->GetNReqStations()) return;
773
774 fhCalDX0->Fill(pTrk->GetFitX(0.));
775 fhCalDY0->Fill(pTrk->GetFitY(0.));
776
777 if (fbBeam && !pTrk->ContainsAddr(CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 5)))
778 return; // request beam counter hit for calibration
779
780 if (fbBeam && fdR0Lim > 0.) // consider only tracks originating from nominal interaction point
781 {
782 fhCalR0->Fill(pTrk->GetR0());
783 if (pTrk->GetR0() > fdR0Lim) return;
784 }
785
786 //if (iOpt > 70) HstDoublets(pTrk); // Fill Doublets histograms
787
788 //double dEvtStart=fEvtHeader->GetEventTime();
789 for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) {
790 CbmTofHit* pHit = pTrk->GetTofHitPointer(iHit);
791 Int_t iDetId = (pHit->GetAddress() & DetMask);
792 Int_t iSmType = CbmTofAddress::GetSmType(iDetId);
793 Int_t iSm = CbmTofAddress::GetSmId(iDetId);
794 Int_t iRpc = CbmTofAddress::GetRpcId(iDetId);
795
796 std::map<UInt_t, UInt_t>::iterator it = fDetIdIndexMap.find(iDetId);
797 if (it == fDetIdIndexMap.end()) continue; // continue for invalid detector index
798 Int_t iDetIndx = it->second; //fDetIdIndexMap[iDetId];
799
800 Int_t iChId = pHit->GetAddress();
801 Int_t iCh = CbmTofAddress::GetChannelId(iChId);
802 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
803 if (NULL == fChannelInfo) {
804 LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
805 continue;
806 }
807 //gGeoManager->FindNode(fChannelInfo->GetX(), fChannelInfo->GetY(), fChannelInfo->GetZ());
808 Double_t hitpos[3];
809 hitpos[0] = pHit->GetX();
810 hitpos[1] = pHit->GetY();
811 hitpos[2] = pHit->GetZ();
812 Double_t hlocal_p[3];
813 //TGeoNode* cNode=
814 //gGeoManager->GetCurrentNode();
815 //gGeoManager->MasterToLocal(hitpos, hlocal_p);
816 fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_p);
817
818 hitpos[0] = pTrk->GetFitX(pHit->GetZ());
819 hitpos[1] = pTrk->GetFitY(pHit->GetZ());
820 Double_t hlocal_f[3];
821 //gGeoManager->MasterToLocal(hitpos, hlocal_f);
822 fTofClusterizer->MasterToLocal(iChId, hitpos, hlocal_f);
823 if (true) { //0==fDigiBdfPar->GetChanOrient(iSmType, iRpc)) { // default
824 fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1]); // transformed into LRF
825 fhCalPos[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1] - hlocal_f[1]); // transformed into LRF
826 }
827 else { // Counters rotated with respect to module axis, strips along x-axis
828 fhCalPosition[iDetIndx]->Fill((Double_t) iCh, hlocal_p[0]); // transformed into LRF
829 fhCalPos[iDetIndx]->Fill((Double_t) iCh, hlocal_p[0] - hlocal_f[0]); // transformed into LRF
830 }
831 //if(fTofFindTracks->GetTtTarg()==0.033) pTrk->SetTt(fTofFindTracks->GetTtTarg()); //for calibration mode 1,2,31
832 //double dTexp = pTrk->GetFitT(pHit->GetZ()); // residuals transformed into LRF
833 //double dTexp = fTrackletTools->GetTexpected(pTrk, pHit->GetAddress()&DetMask, pHit, fTofFindTracks->GetTtLight());
834 double dTexp = fTrackletTools->GetTexpected(pTrk, -1, pHit, fTofFindTracks->GetTtLight());
835 fhCalTOff[iDetIndx]->Fill((Double_t) iCh, pHit->GetTime() - dTexp);
836 fhCalXYTOff[iDetIndx]->Fill((Double_t) iCh, hlocal_p[1],
837 pHit->GetTime() - dTexp); // residuals transformed into LRF
838 //fhCalTOff[iDetIndx]->Fill((Double_t)iCh,fTrackletTools->GetTdif(pTrk, iDetId, pHit));
839
840 Int_t iEA = pTrk->GetTofHitIndex(iHit);
841 Int_t iTSA = fTofFindTracks->GetTofHitIndex(iEA);
842
843 if (iTSA > fTofDigiMatchColl->GetEntriesFast()) {
844 LOG(error) << " Inconsistent DigiMatches for Hitind " << iTSA
845 << ", TClonesArraySize: " << fTofDigiMatchColl->GetEntriesFast();
846 }
847 CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(iTSA);
848
849 Double_t hlocal_d[3];
850 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink += 2) { // loop over digis
851 CbmLink L0 = digiMatch->GetLink(iLink);
852 Int_t iDigInd0 = L0.GetIndex();
853 Int_t iDigInd1 = (digiMatch->GetLink(iLink + 1)).GetIndex();
854
855 const CbmTofDigi* tDigi0 = NULL;
856 const CbmTofDigi* tDigi1 = NULL;
857 if (tEvent != NULL) { //disable
858 LOG(debug) << "Locate MatchDigiInd " << iDigInd0 << " and " << iDigInd1 << " in CalDigiVec of size "
859 << fTofCalDigiVec->size();
860 // <<" in current event not implemented";
861 //continue;
862 tDigi0 = &(fTofCalDigiVec->at(iDigInd0));
863 tDigi1 = &(fTofCalDigiVec->at(iDigInd1));
864 }
865 else { // event wise entries (TS mode)
866 LOG(debug) << "TS mode: locate MatchDigiInd " << iDigInd0 << " and " << iDigInd1 << " in CalDigiVec of size "
867 << fTofCalDigiVec->size();
868 if (fair::Logger::Logging(fair::Severity::debug2)) {
869 for (UInt_t iDigi = 0; iDigi < fTofCalDigiVec->size(); iDigi++) {
870 tDigi0 = &(fTofCalDigiVec->at(iDigi));
871 LOG(debug) << Form("#%u: TSRCS %d%d%d%2d%d", iDigi, (Int_t) tDigi0->GetType(), (Int_t) tDigi0->GetSm(),
872 (Int_t) tDigi0->GetRpc(), (Int_t) tDigi0->GetChannel(), (Int_t) tDigi0->GetSide());
873 }
874 }
875 tDigi0 = &(fTofCalDigiVec->at(iDigInd0));
876 tDigi1 = &(fTofCalDigiVec->at(iDigInd1));
877 //tDigi0 = fDigiMan->Get<CbmTofDigi>(iDigInd0);
878 //tDigi1 = fDigiMan->Get<CbmTofDigi>(iDigInd1);
879 }
880
881 Int_t iCh0 = tDigi0->GetChannel();
882 Int_t iSide0 = tDigi0->GetSide();
883
884 LOG(debug) << "Fill Walk for Hit Ind " << iEA << ", " << iTSA
885 << Form(", TSRC %d%d%d%2d, DigiInd %2d, %2d", iSmType, iSm, iRpc, iCh, iDigInd0, iDigInd1)
886 << Form(", TSRCS %d%d%d%2d%d %d%d%d%2d%d", (Int_t) tDigi0->GetType(), (Int_t) tDigi0->GetSm(),
887 (Int_t) tDigi0->GetRpc(), (Int_t) tDigi0->GetChannel(), (Int_t) tDigi0->GetSide(),
888 (Int_t) tDigi1->GetType(), (Int_t) tDigi1->GetSm(), (Int_t) tDigi1->GetRpc(),
889 (Int_t) tDigi1->GetChannel(), (Int_t) tDigi1->GetSide());
890
891 if (iDetIndx > (Int_t) fhCalWalk.size()) {
892 LOG(error) << "Invalid DetIndx " << iDetIndx;
893 continue;
894 }
895 if (iCh0 > (Int_t) fhCalWalk[iDetIndx].size()) {
896 LOG(error) << "Invalid Ch0 " << iCh0 << " for detIndx" << iDetIndx << ", size " << fhCalWalk[iDetIndx].size();
897 continue;
898 }
899 if (iSide0 > (Int_t) fhCalWalk[iDetIndx][iCh0].size()) {
900 LOG(error) << "Invalid Side0 " << iSide0 << " for "
901 << " for detIndx" << iDetIndx << ", size " << fhCalWalk[iDetIndx][iCh0].size();
902 continue;
903 }
904
905 Int_t iCh1 = tDigi1->GetChannel();
906 Int_t iSide1 = tDigi1->GetSide();
907 if (iCh1 > (Int_t) fhCalWalk[iDetIndx].size()) {
908 LOG(error) << "Invalid Ch1 " << iCh1 << " for "
909 << " for detIndx" << iDetIndx;
910 continue;
911 }
912 if (iSide1 > (Int_t) fhCalWalk[iDetIndx][iCh1].size()) {
913 LOG(error) << "Invalid Side1 " << iSide1 << " for "
914 << " for detIndx" << iDetIndx;
915 continue;
916 }
917
918 if (iCh0 != iCh1 || iSide0 == iSide1) {
919 LOG(fatal) << "Invalid digi pair for TSR " << iSmType << iSm << iRpc
920 << Form(" Ch %2d %2d side %d %d", iCh0, iCh1, iSide0, iSide1) << " in event " << NevtT;
921 continue;
922 }
923
924 hlocal_d[1] =
925 -0.5 * ((1. - 2. * tDigi0->GetSide()) * tDigi0->GetTime() + (1. - 2. * tDigi1->GetSide()) * tDigi1->GetTime())
926 * fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
927
928 if (TMath::Abs(hlocal_d[1] - hlocal_p[1]) > 10.) {
929 LOG(warn) << "CMPY for TSRC " << iSmType << iSm << iRpc << iCh0 << ": " << hlocal_f[1] << ", " << hlocal_p[1]
930 << ", " << hlocal_d[1] << ", TOT: " << tDigi0->GetTot() << " " << tDigi1->GetTot();
931 }
932 Int_t iWalkMode = 0; //(iOpt - iOpt % 100) / 100;
933 switch (iWalkMode) {
934 case 1:
935 fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(
936 tDigi0->GetTot(),
937 tDigi0->GetTime() + (1. - 2. * tDigi0->GetSide()) * hlocal_d[1] / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)
938 - pTrk->GetFitT(pHit->GetZ()) //-fTrackletTools->GetTexpected(pTrk, iDetId, pHit)
939 + fTofFindTracks->GetTOff(iDetId)
940 + 2. * (1. - 2. * tDigi0->GetSide()) * (hlocal_d[1] - hlocal_f[1])
941 / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
942 /*
943 LOG(info)<<"TSRCS "<<iSmType<<iSm<<iRpc<<iCh<<iSide0<<Form(": digi0 %f, ex %f, prop %f, Off %f, res %f",
944 tDigi0->GetTime(),
945 fTrackletTools->GetTexpected(pTrk, iDetId, pHit) ,
946 fTofFindTracks->GetTOff(iDetId),
947 (1.-2.*tDigi0->GetSide())*hlocal_f[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc),
948 tDigi0->GetTime()-fTrackletTools->GetTexpected(pTrk, iDetId, pHit)
949 -(1.-2.*tDigi0->GetSide())*hlocal_f[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc));
950 */
951
952 fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(
953 tDigi1->GetTot(),
954 tDigi1->GetTime() + (1. - 2. * tDigi1->GetSide()) * hlocal_d[1] / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc)
955 - pTrk->GetFitT(pHit->GetZ()) //-fTrackletTools->GetTexpected(pTrk, iDetId, pHit)
956 + fTofFindTracks->GetTOff(iDetId)
957 + 2. * (1. - 2. * tDigi1->GetSide()) * (hlocal_d[1] - hlocal_f[1])
958 / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
959 break;
960
961 case 0: {
962 Double_t dDeltaT =
963 0.5 * (tDigi0->GetTime() + tDigi1->GetTime()) + fTofFindTracks->GetTOff(iDetId) - pHit->GetTime();
964
965 fhCalWalkAv[iDetIndx]->Fill(tDigi0->GetTot(), dDeltaT);
966 fhCalWalkAv[iDetIndx]->Fill(tDigi1->GetTot(), dDeltaT);
967 fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi0->GetTot(), dDeltaT);
968 fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi1->GetTot(), dDeltaT);
969 if (iSmType == 5 || iSmType == 8) { // symmetrize for Pad counters
970 fhCalWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi1->GetTot(), dDeltaT);
971 fhCalWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi0->GetTot(), dDeltaT);
972 }
973 Double_t dDeltaDT = (hlocal_d[1] - hlocal_p[1]) / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
974 fhCalDtWalk[iDetIndx][iCh0][iSide0]->Fill(tDigi0->GetTot(), (2. * tDigi0->GetSide() - 1) * dDeltaDT);
975 fhCalDtWalk[iDetIndx][iCh1][iSide1]->Fill(tDigi1->GetTot(), (2. * tDigi1->GetSide() - 1) * dDeltaDT);
976 } break;
977 }
978 }
979 }
980}
981
982Bool_t CbmTofCalibrator::UpdateCalHist(Int_t iOpt)
983{
984 // get current calibration histos
985 TString CalFileName = fTofClusterizer->GetCalParFileName();
986 LOG(info) << "CbmTofCalibrator:: update histos from "
987 << "file " << CalFileName << " with option " << iOpt;
988 int iOpt0 = iOpt % 10;
989 int iOpt1 = (iOpt - iOpt0) / 10;
990
992 TFile* oldFile = gFile;
993 TDirectory* oldDir = gDirectory;
994
995 TFile* fCalParFile = new TFile(CalFileName, "update");
996 if (NULL == fCalParFile) {
997 LOG(warn) << "Could not open TofClusterizer calibration file " << CalFileName;
998 if (iOpt0 == 9) { //modify reference file name
999 int iCalMode = CalFileName.Index("tofClust") - 3;
1000 CalFileName(iCalMode) = '3';
1001 LOG(info) << "Modified CalFileName = " << CalFileName;
1002 fCalParFile = new TFile(CalFileName, "update");
1003 if (NULL == fCalParFile) LOG(fatal) << "Could not open TofClusterizer calibration file " << CalFileName;
1004 }
1005 }
1006 assert(fCalParFile);
1007 ReadHist(fCalParFile);
1008 if (kTRUE) { // all calibration modes
1009 const Double_t MINCTS = 100.; //FIXME, numerical constant in code
1010 // modify calibration histograms
1011 // check for beam counter
1012 Double_t dBeamTOff = 0.;
1013 for (Int_t iDetIndx = 0; iDetIndx < fDigiBdfPar->GetNbDet(); iDetIndx++) {
1014 Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
1015 Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
1016 if (5 == iSmType) {
1017 if (fhCalTOff[iDetIndx]->GetEntries() > MINCTS) {
1018 TH1* hBy = (TH1*) fhCalTOff[iDetIndx]->ProjectionY();
1019 // Fit gaussian around peak value
1020 Double_t dFMean = hBy->GetBinCenter(hBy->GetMaximumBin());
1021 Double_t dFLim = 0.5; // CAUTION, fixed numeric value
1022 Double_t dBinSize = hBy->GetBinWidth(1);
1023 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1024 TFitResultPtr fRes = hBy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1025 dBeamTOff = fRes->Parameter(1); //overwrite mean
1026 LOG(info) << "Found beam counter with average TOff = " << dBeamTOff;
1027 }
1028 else {
1029 LOG(info) << "Beam counter has too few entries: " << fhCalTOff[iDetIndx]->GetEntries();
1030 }
1031 break;
1032 }
1033 }
1034 if (dBeamTOff == 0.) LOG(warn) << "No beam counter found";
1035
1036 for (Int_t iDetIndx = 0; iDetIndx < fDigiBdfPar->GetNbDet(); iDetIndx++) {
1037 Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
1038 // Int_t iSmAddr = iUniqueId & DetMask;
1039 Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
1040 Int_t iSm = CbmTofAddress::GetSmId(iUniqueId);
1041 Int_t iRpc = CbmTofAddress::GetRpcId(iUniqueId);
1042 if (NULL == fhCorTOff[iDetIndx]) {
1043 LOG(warn) << "hCorTOff for TSR " << iSmType << iSm << iRpc << " not available";
1044 continue;
1045 }
1046
1047 switch (iOpt0) {
1048 case 0: // none
1049 {
1050 double* dRes;
1051 switch (iOpt1) {
1052 case 10: {
1053 LOG(info) << "Equalize velocities with const threshold counterwise for " << fhCorTOff[iDetIndx]->GetName()
1054 << ", MeanTOff " << fhCorTOff[iDetIndx]->GetMean();
1055 TString hname2 = Form("cal_SmT%d_sm%03d_rpc%03d_TOff", iSmType, iSm, iRpc);
1056 //TH2* h2 = (TH2*) gROOT->FindObjectAny(hname2);
1057 TH2* h2 = fhCalTOff[iDetIndx];
1058 TH2* h2W = fhCalTofOff[iDetIndx];
1059 if (NULL == h2) {
1060 LOG(warn) << "Calibration data not available from " << hname2;
1061 continue;
1062 }
1063
1064 double dResIni[1] = {0.};
1065 dRes = &dResIni[0];
1066 Int_t NStrips = h2->GetNbinsX();
1067 TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
1068 TH1* h2y = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
1069 double dEdge = 0.;
1070 if (iSmType != 5) {
1071 dRes = find_tofedge((const char*) (h2Wy->GetName()), fTofClusterizer->GetEdgeThr(),
1072 fTofClusterizer->GetEdgeLen());
1073 dEdge = dRes[0];
1074 }
1075 else {
1076 dEdge = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
1077 }
1078
1079 /*
1080 LOG(info) << h2Wy->GetName() << ": Coarse Max at " << h2Wy->GetMaximumBin()
1081 << ", " << dEdge << ", " << h2y->GetBinCenter(1)
1082 << ", " << h2y->GetBinCenter(h2y->GetNbinsX());
1083 */
1084 if (dEdge < h2y->GetBinCenter(1) || dEdge > h2y->GetBinCenter(h2y->GetNbinsX())) {
1085 dRes[0] = dEdge; // - fTofClusterizer->GetEdgeTbias();
1086 LOG(info) << h2Wy->GetName() << ": Coarse TOff shift by " << dRes[0];
1087 }
1088 else {
1089 if (iSmType != 5) {
1090 dRes = find_tofedge((const char*) (h2y->GetName()), fTofClusterizer->GetEdgeThr(),
1091 fTofClusterizer->GetEdgeLen());
1092 if (TMath::Abs(dRes[0]) < 0.5) {
1093 //LOG(info) << "Now fit edge above noise, dT " << dRes[0] <<", T "<< iSmType;
1094 double dTfind = dRes[0];
1095 double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
1096 dRes = fit_tofedge((const char*) (h2y->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
1097 if (std::isnan(dRes[0])) dRes[0] = dTfind;
1098 }
1099 }
1100 }
1101 //LOG(info) << h2y->GetName() << ": shift TOff by " << dRes[0];
1102 fhCalChannelDt->Fill(dRes[0]);
1103 for (Int_t i = 0; i < NStrips; i++) {
1104 fhCorTOff[iDetIndx]->SetBinContent(i + 1, fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0]);
1105 }
1106 } break;
1107
1108 case 15:
1109 case 12:
1110 case 11: {
1111 LOG(info) << "Equalize velocities with chi2 deviations for " << fhCorTOff[iDetIndx]->GetName();
1112 TString hname2 = Form("cal_SmT%d_sm%03d_rpc%03d_TOff", iSmType, iSm, iRpc);
1113
1114 double dTRefSum = 0.;
1115 double dRefNormSum = 0.;
1116 double dNorm = 0.;
1117
1118 TH2* h2W = fhCalTofOff[iDetIndx];
1119
1120 if (h2W->GetEntries() < MINCTS) continue;
1121
1122 TH2* h2 = fhCalTOff[iDetIndx];
1123 Int_t NStrips = h2->GetNbinsX();
1124 TH1* h2yAv;
1125 double dResIni[1] = {0.};
1126 dRes = &dResIni[0];
1127
1128 TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
1129 h2yAv = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
1130 double dEdge = 0.;
1131
1132 double dMean = TruncatedMeanY(h2W); //h2Wy->GetMean();
1133 if (iSmType == 5) {
1134 dEdge = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
1135 dNorm = h2Wy->GetBinContent(h2Wy->GetMaximumBin());
1136 }
1137 else {
1138 dRes = find_tofedge((const char*) (h2Wy->GetName()), fTofClusterizer->GetEdgeThr(),
1139 fTofClusterizer->GetEdgeLen());
1140 dEdge = dRes[0];
1141 if ((dMean - dEdge) > 10.) { //FIXME: constant in code
1142 LOG(warn) << h2Wy->GetName() << ": wrong edge detected " << dEdge << ", Mean " << dMean;
1143 dEdge = dMean;
1144 }
1145 }
1146 /*
1147 LOG(info) << h2Wy->GetName() << ": Coarse shift at " << h2Wy->GetMaximumBin()
1148 << ", Edge " << dEdge
1149 << ", " << h2yAv->GetBinCenter(1)
1150 << ", " << h2yAv->GetBinCenter(h2yAv->GetNbinsX());
1151 */
1152 if (dEdge < h2yAv->GetBinCenter(1) || dEdge > h2yAv->GetBinCenter(h2yAv->GetNbinsX())) {
1153 dRes[0] = dEdge; // - fTofClusterizer->GetEdgeTbias();
1154 //LOG(info) << h2Wy->GetName() << ": Coarse TOff shift by " << dRes[0];
1155 }
1156 else {
1157 if (iSmType != 5) {
1158 dRes = find_tofedge((const char*) (h2yAv->GetName()), fTofClusterizer->GetEdgeThr(),
1159 fTofClusterizer->GetEdgeLen());
1160 if (TMath::Abs(dRes[0]) < dValidEdge) {
1161 double dTfind = dRes[0];
1162 double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
1163 dRes = fit_tofedge((const char*) (h2yAv->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
1164 if (std::isnan(dRes[0])) dRes[0] = dTfind;
1165 LOG(info) << h2yAv->GetName() << ": fitted AvEdge up to " << dTmax << ", Res " << dRes[0]
1166 << ", FindE " << dTfind;
1167 }
1168 }
1169 }
1170
1171 if (iSmType != 5) {
1172 fhCalCounterDt->Fill(dRes[0] - fTofClusterizer->GetEdgeTbias());
1173 if (TMath::Abs(dRes[0] - fTofClusterizer->GetEdgeTbias()) > -fhCalCounterDt->GetBinCenter(1)) {
1174 LOG(warn) << "Channel ooCR: " << fhCorTOff[iDetIndx]->GetName() << ", " << dRes[0] << ", "
1175 << fTofClusterizer->GetEdgeTbias() << ", " << -fhCalCounterDt->GetBinCenter(1);
1176 }
1177 }
1178 else {
1179 fhCalCounterDt->Fill(dEdge);
1180 }
1181
1182 const double dResAv = dRes[0];
1183 const int nValues = MaxShift * 2 + 1;
1184 const double dRangeFac = 1.;
1185 double chi2Shift[NStrips][nValues];
1186 double chi2Minimum[NStrips];
1187 double chi2MinShift[NStrips];
1188 double chi2LowMinNeighbor[NStrips];
1189 double chi2LowMinNeighborShift[NStrips];
1190 for (Int_t i = 0; i < NStrips; i++) {
1191 double dTminShift = 0.;
1192 double dEdgei = 0.;
1193 double dMeani = 0.;
1194 TH1* h2y = h2->ProjectionY(Form("%s_py%d", h2->GetName(), i), i + 1, i + 1);
1195 if (h2y->GetEntries() < MINCTS) continue;
1196 LOG(info) << h2y->GetName() << ": inspect with iOpt1 = " << iOpt1;
1197 if (kTRUE) {
1198 if (iOpt1 == 11 || iOpt1 == 15) {
1199 //check peak position
1200 TH1* h2Wyi = h2W->ProjectionY(Form("%s_py%d", h2W->GetName(), i), i + 1, i + 1);
1201 if (iSmType == 5) {
1202 dEdgei = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
1203 dNorm = h2Wy->GetBinContent(h2Wy->GetMaximumBin());
1204 }
1205 else {
1206 dRes = find_tofedge((const char*) (h2Wyi->GetName()), fTofClusterizer->GetEdgeThr(),
1207 fTofClusterizer->GetEdgeLen());
1208 dEdgei = dRes[0];
1209 }
1210
1211 if (dEdgei < dRangeFac * h2y->GetBinCenter(1)
1212 || dEdgei > dRangeFac * h2y->GetBinCenter(h2y->GetNbinsX())) {
1213 dMeani = h2Wyi->GetMean();
1214 LOG(info) << h2Wyi->GetName() << " ooR: " << dRes[0] << ", E " << dEdgei << ", Av " << dResAv
1215 << ", Meani " << dMeani << ", MeanY " << dMean;
1216 if ((dMeani - dEdgei) > 10.) { //FIXME: constant in code
1217 LOG(warn) << h2Wyi->GetName() << ": invalid edge detected " << dEdgei << ", Mean " << dMeani;
1218 dRes[0] = dMeani - dMean;
1219 }
1220 //dRes[0]=dEdgei
1221 }
1222 else {
1223 if (iSmType != 5) {
1224 dRes = find_tofedge((const char*) (h2y->GetName()), fTofClusterizer->GetEdgeThr(),
1225 fTofClusterizer->GetEdgeLen());
1226 LOG(info) << h2y->GetName() << " CheckValidEdge: " << dRes[0] << ", " << dValidEdge;
1227 if (TMath::Abs(dRes[0]) < dValidEdge) {
1228 switch (iOpt1) {
1229 case 11: {
1230 double dTfind = dRes[0];
1231 double dTmax = dTfind + fTofClusterizer->GetEdgeFrange();
1232 dRes = fit_tofedge((const char*) (h2y->GetName()), dTmax, fTofClusterizer->GetEdgeThr());
1233 if (std::isnan(dRes[0])) dRes[0] = dTfind;
1234 LOG(info) << h2y->GetName() << ": fit edge up to " << dTmax << ", Res " << dRes[0]
1235 << ", FindE " << dTfind;
1236 } break;
1237 case 15:
1238 chi2Minimum[i] = 1.E6;
1239 chi2LowMinNeighbor[i] = 1.E6;
1240 int idx = 0;
1241 for (int j = -MaxShift; j < MaxShift + 1; j++) {
1242 chi2Shift[i][idx] = CalcChi2(h2yAv, h2y, j);
1243 if (chi2Shift[i][idx] < chi2Minimum[i]) {
1244 chi2LowMinNeighbor[i] = chi2Minimum[i];
1245 chi2LowMinNeighborShift[i] = chi2MinShift[i];
1246 chi2Minimum[i] = chi2Shift[i][idx];
1247 chi2MinShift[i] = j * h2y->GetBinWidth(1);
1248 }
1249 else {
1250 if (chi2Shift[i][idx] < chi2LowMinNeighbor[i]) {
1251 chi2LowMinNeighbor[i] = chi2Shift[i][idx];
1252 chi2LowMinNeighborShift[i] = j * h2y->GetBinWidth(1);
1253 }
1254 }
1255 idx++;
1256 }
1257 double chi2sum = chi2Minimum[i] + chi2LowMinNeighbor[i];
1258 dTminShift = (chi2MinShift[i] * (chi2sum - chi2Minimum[i])
1259 + chi2LowMinNeighborShift[i] * (chi2sum - chi2LowMinNeighbor[i]))
1260 / chi2sum;
1261 //dTminShift=chi2MinShift[i]; // take Min only, for debugging
1262 LOG(info) << h2y->GetName() << ": Chi2 shift " << dTminShift << ", Res " << dRes[0]
1263 << ", ResAv " << dResAv;
1264 dRes[0] = dResAv;
1265 break;
1266 }
1267 }
1268 else {
1269 LOG(info) << h2y->GetName() << ": stick to coarse offset " << dEdgei << ", Res " << dRes[0];
1270 dRes[0] = dEdgei; // use offset from wide histo
1271 }
1272 }
1273 else { // beam counter
1274 dEdgei = h2y->GetBinCenter(h2y->GetMaximumBin());
1275 dNorm = h2y->GetBinContent(h2y->GetMaximumBin());
1276 }
1277
1278 //LOG(info) << h2Wyi->GetName() << " TminShift: " << dTminShift;
1279 }
1280 }
1281 }
1282 if (iSmType == 5) {
1283 LOG(info) << "Update StartCounterCalib " << i << ": " << dRes[0] << ", " << dBeamTOff << ", "
1284 << dTminShift << ", " << dEdgei << ", N " << dNorm;
1285 dRes[0] = -dEdgei;
1286 }
1287 else {
1288 dTminShift -= fTofClusterizer->GetEdgeTbias(); // shift to physical scale
1289 }
1290 fhCalChannelDt->Fill(dRes[0] + dTminShift);
1291 if (TMath::Abs(dRes[0] + dTminShift) > -fhCalChannelDt->GetBinCenter(1)) {
1292 LOG(warn) << "Channel ooHR: " << h2Wy->GetName() << i << ", " << dRes[0] << ", " << dTminShift;
1293 }
1294 double dVal = fhCorTOff[iDetIndx]->GetBinContent(i + 1);
1295 if (std::isnan(dVal)) fhCorTOff[iDetIndx]->SetBinContent(i + 1, 0.); // fix NANs
1296
1297 LOG(info) << h2Wy->GetName() << i << " CorTOff: " << dRes[0] << ", " << dResAv << ", " << dTminShift
1298 << ", " << fhCorTOff[iDetIndx]->GetBinContent(i + 1) << " -> "
1299 << fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0] + dTminShift;
1300
1301 fhCorTOff[iDetIndx]->SetBinContent(i + 1,
1302 (fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dRes[0] + dTminShift));
1303 if (iSmType == 5) {
1304 dTRefSum += fhCorTOff[iDetIndx]->GetBinContent(i + 1) * dNorm;
1305 dRefNormSum += dNorm;
1306 }
1307 //LOG(info)<<fhCorTOff[iDetIndx]->GetName() << ": new TOff " << fhCorTOff[iDetIndx]->GetBinContent(i + 1);
1308 } //for (Int_t i = 0; i < NStrips; i++)
1309 if (iSmType == 5) { // avoid runaway
1310 double dTRefAv = dTRefSum / dRefNormSum;
1311 LOG(info) << "Shift average ref time to 0: " << dTRefAv;
1312 for (Int_t i = 0; i < NStrips; i++) {
1313 fhCorTOff[iDetIndx]->SetBinContent(i + 1, (fhCorTOff[iDetIndx]->GetBinContent(i + 1) - dTRefAv));
1314 }
1315 }
1316 } break;
1317
1318 case 14:
1319 case 13: {
1320 LOG(info) << "TofCalibrator: skip TOff, calibrate hit positions, check YBox, Opt = " << iOpt;
1321 } break;
1322 case 16: {
1323 LOG(info) << "TofCalibrator: determine counter time offsets with cosmics, Opt = " << iOpt;
1324 TH2* h2 = fhCalTOff[iDetIndx];
1325 TH2* h2W = fhCalTofOff[iDetIndx];
1326 if (NULL == h2) {
1327 LOG(warn) << "Calibration data not available for detx " << iDetIndx;
1328 continue;
1329 }
1330 Int_t NStrips = h2->GetNbinsX();
1331 TH1* h2Wy = h2W->ProjectionY(Form("%s_py", h2W->GetName()), 1, NStrips);
1332 TH1* h2y = h2->ProjectionY(Form("%s_py", h2->GetName()), 1, NStrips);
1333
1334 Double_t dCtsWy = h2Wy->GetEntries();
1335 Double_t dDt = 0;
1336 if (dCtsWy > MINCTS) {
1337 double dFMean = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
1338 double dFLim = 0.5; // CAUTION, fixed numeric value
1339 double dBinSize = h2Wy->GetBinWidth(1);
1340 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1341 TFitResultPtr fRes = h2Wy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1342 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1343 dDt = fRes->Parameter(1);
1344 if (TMath::Abs(dDt) < 0.9 * h2y->GetXaxis()->GetXmax()) {
1345 dFMean = h2y->GetBinCenter(h2y->GetMaximumBin());
1346 dFLim = 0.5; // CAUTION, fixed numeric value
1347 dBinSize = h2y->GetBinWidth(1);
1348 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1349 fRes = h2y->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1350 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1351 dDt = fRes->Parameter(1);
1352 }
1353 }
1354 fhCalCounterDt->Fill(dDt);
1355 }
1356 LOG(info) << "Update cor hist " << fhCorTOff[iDetIndx]->GetName() << " by " << dDt << ", " << dBeamTOff;
1357 for (int iCh = 0; iCh < fhCorTOff[iDetIndx]->GetNbinsX(); iCh++) {
1358 fhCalChannelDt->Fill(dDt);
1359 fhCorTOff[iDetIndx]->SetBinContent(iCh + 1, fhCorTOff[iDetIndx]->GetBinContent(iCh + 1) + dDt);
1360 }
1361 }
1362 } break;
1363
1364 case 17: {
1365 LOG(info) << "TofCalibrator: determine channel time offsets with cosmics, Opt = " << iOpt;
1366 TH2* h2 = fhCalTOff[iDetIndx];
1367 TH2* h2W = fhCalTofOff[iDetIndx];
1368 if (NULL == h2) {
1369 LOG(warn) << "Calibration data not available for detx " << iDetIndx;
1370 continue;
1371 }
1372 Int_t NStrips = h2->GetNbinsX();
1373 for (int i = 0; i < NStrips; i++) {
1374 TH1* h2Wy = h2W->ProjectionY(Form("%s_py_%02d", h2W->GetName(), i), i + 1, i + 1);
1375 TH1* h2y = h2->ProjectionY(Form("%s_py_%02d", h2->GetName(), i), i + 1, i + 1);
1376
1377 Double_t dCtsWy = h2Wy->GetEntries();
1378 Double_t dDt = 0;
1379 if (dCtsWy > MINCTS) {
1380 double dFMean = h2Wy->GetBinCenter(h2Wy->GetMaximumBin());
1381 double dFLim = 0.5; // CAUTION, fixed numeric value
1382 double dBinSize = h2Wy->GetBinWidth(1);
1383 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1384 TFitResultPtr fRes = h2Wy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1385 Int_t iFitStatus = fRes;
1386 //LOG(warn)<<h2Wy->GetName() << ": "<< iFitStatus << " " <<fRes;
1387 if (iFitStatus != -1)
1388 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1389 dDt = fRes->Parameter(1);
1390 //LOG(warn)<<h2Wy->GetName() << " Dt "<< dDt;
1391 if (h2y->GetEntries() > MINCTS)
1392 if (TMath::Abs(dDt) < 0.9 * h2y->GetXaxis()->GetXmax()) {
1393 dFMean = h2y->GetBinCenter(h2y->GetMaximumBin());
1394 dFLim = 0.5; // CAUTION, fixed numeric value
1395 dBinSize = h2y->GetBinWidth(1);
1396 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1397 fRes = h2y->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1398 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1399 dDt = fRes->Parameter(1);
1400 }
1401 }
1402 fhCalChannelDt->Fill(dDt);
1403 fhCorTOff[iDetIndx]->SetBinContent(i + 1, fhCorTOff[iDetIndx]->GetBinContent(i + 1) + dDt);
1404 }
1405 }
1406 else {
1407 LOG(warn) << "Too few counts in " << h2Wy->GetName() << ": " << dCtsWy;
1408 }
1409 }
1410 } break;
1411
1412 case 18: {
1413 LOG(info) << "TofCalibrator: determine channel walks with cosmics, Opt = " << iOpt << " for TSR "
1414 << iSmType << iSm << iRpc;
1415 //if (iSmType == 5) continue; // no walk correction for beam counter up to now
1416 const Double_t MinCounts = 10.;
1417 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1418 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
1419 TH1* hCY = fhCalPos[iDetIndx]->ProjectionY(Form("%s_py_%d", fhCalPos[iDetIndx]->GetName(), iCh),
1420 iCh + 1, iCh + 1);
1421 if (hCY->GetEntries() > 1)
1422 fhCalChannelDy->Fill(hCY->GetRMS() / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc));
1423 //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
1424 TProfile* hpW0 = fhCalWalk[iDetIndx][iCh][0]->ProfileX(); // mean deviation
1425 TH1* hCW0 = fhCalWalk[iDetIndx][iCh][0]->ProjectionX(); // contributing counts
1426 TProfile* hpW1 = fhCalWalk[iDetIndx][iCh][1]->ProfileX(); // mean deviation
1427 TH1* hCW1 = fhCalWalk[iDetIndx][iCh][1]->ProjectionX(); // contributing counts
1428 if (hCW0->GetEntries() == 0 || hCW1->GetEntries() == 0) {
1429 //LOG(info) << "No entries in " << hCW0->GetName();
1430 continue;
1431 }
1432 Double_t dCorT = 0;
1433 for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][0]->GetNbinsX(); iBin++) {
1434 Double_t dCts0 = hCW0->GetBinContent(iBin + 1);
1435 Double_t dCts1 = hCW1->GetBinContent(iBin + 1);
1436 Double_t dWOff0 = fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin + 1); // current value
1437 Double_t dWOff1 = fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin + 1); // current value
1438 if (iBin > 0 && dCts0 == 0 && dCts1 == 0) {
1439 fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1,
1440 fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin));
1441 fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1,
1442 fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin));
1443 }
1444 dCorT = 0.;
1445 if (dCts0 > MinCounts && dCts1 > MinCounts) {
1446 dCorT = 0.5 * (hpW0->GetBinContent(iBin + 1) + hpW1->GetBinContent(iBin + 1));
1447 fhCalChannelDt->Fill(dCorT);
1448 }
1449 fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1, dWOff0 + dCorT); //set new value
1450 fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1, dWOff1 + dCorT); //set new value
1451 }
1452 }
1453 } break;
1454
1455 default:
1456 LOG(info) << "TofCalibrator: unknown option " << iOpt;
1457 return kTRUE;
1458 break;
1459 } // switch iOpt1 end
1460
1461 // re-adjust positions
1462 if (fhCalPosition[iDetIndx]->GetEntries() == 0) continue;
1463 double dYShift = 0;
1464 TH1* hCalP = fhCalPosition[iDetIndx]->ProjectionX();
1465 TH1* hCalPos_py = fhCalPosition[iDetIndx]->ProjectionY(); // full counter
1466 dYShift = hCalPos_py->GetMean();
1467 int iChId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
1468 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
1469 if (NULL == fChannelInfo) LOG(fatal) << Form("invalid ChannelInfo for 0x%08x", iChId);
1470 double dLen = fChannelInfo->GetSizey() * fTofClusterizer->GetModifySigvel();
1471 fTofClusterizer->fit_ybox(hCalPos_py, dLen);
1472 TF1* ff = hCalPos_py->GetFunction("YBox"); // box fit for monitoring signal velocity
1473 if (NULL != ff) {
1474 LOG(info) << "FRes YBox " << hCalPos_py->GetEntries() << " entries in TSR " << iSmType << iSm << iRpc
1475 << ", chi2 " << ff->GetChisquare() / ff->GetNDF()
1476 << Form(", striplen (%5.2f): %7.2f +/- %5.2f, pos "
1477 "res %5.2f +/- %5.2f at y_cen = %5.2f +/- %5.2f",
1478 dLen * fTofClusterizer->GetModifySigvel(), 2. * ff->GetParameter(1),
1479 2. * ff->GetParError(1), ff->GetParameter(2), ff->GetParError(2), ff->GetParameter(3),
1480 ff->GetParError(3));
1481 //dYShift = ff->GetParameter(3); // update common shift by box fit
1482 if (iOpt1 == 14) { // Update signal velocity
1483 if (hSvel[iSmType] != NULL) {
1484 int iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
1485 double dSscal = hSvel[iSmType]->GetBinContent(iSm * iNbRpc + iRpc + 1);
1486 double dRatio = dLen * fTofClusterizer->GetModifySigvel() / 2. / ff->GetParameter(1);
1487 LOG(info) << "Modify SigVel modifier for TSR " << iSmType << iSm << iRpc << ": " << dSscal << ", "
1488 << dRatio << " -> " << dSscal * dRatio;
1489 dSscal = dSscal * dRatio;
1490 if (dSscal < 0.8) dSscal = 0.8;
1491 if (dSscal > 1.2) dSscal = 1.2;
1492 hSvel[iSmType]->SetBinContent(iSm * iNbRpc + iRpc + 1,
1493 dSscal); // does not work in input file directory
1494 }
1495 }
1496 }
1497
1498 double dThr = hCalPos_py->GetBinContent(hCalPos_py->GetMaximumBin()) / 10;
1499 if (dThr < 3.) dThr = 3.;
1500 dRes = fTofClusterizer->find_yedges((const char*) (hCalPos_py->GetName()), dThr, fChannelInfo->GetSizey());
1501 if (TMath::Abs(dRes[0] - fChannelInfo->GetSizey()) / fChannelInfo->GetSizey() < 0.1) {
1502 dYShift = dRes[1];
1503 }
1504
1505 if (iOpt1 != 16 && iOpt1 != 12) {
1506 for (Int_t iBin = 0; iBin < fhCorPos[iDetIndx]->GetNbinsX(); iBin++) {
1507 Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1508 Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
1509 if (dCtsP > MINCTS) {
1510 TH1* hPos_py = fhCalPosition[iDetIndx]->ProjectionY(
1511 Form("%s_py_%d", fhCalPosition[iDetIndx]->GetName(), iBin), iBin + 1, iBin + 1);
1512 dYShift = hPos_py->GetMean();
1513 dThr = hPos_py->GetBinContent(hPos_py->GetMaximumBin()) / 10;
1514 if (dThr < 3.) dThr = 3.;
1515 dRes = fTofClusterizer->find_yedges((const char*) (hPos_py->GetName()), dThr, fChannelInfo->GetSizey());
1516
1517 LOG(info) << Form("EdgeY for %s, TSR %d%d%d: DY %5.2f, Len %5.2f, Size %5.2f ", hPos_py->GetName(),
1518 iSmType, iSm, iRpc, dRes[1], dRes[0],
1519 fChannelInfo->GetSizey() * fTofClusterizer->GetModifySigvel());
1520
1521 if (TMath::Abs(dRes[0] - fChannelInfo->GetSizey()) / fChannelInfo->GetSizey() < 0.1) {
1522 dYShift = dRes[1];
1523 }
1524 else {
1525 dYShift = hPos_py->GetMean();
1526 }
1527 // apply correction
1528 //LOG(info)<<Form("UpdateCalPos TSRC %d%d%d%2d: %6.2f -> %6.2f ",iSmType,iSm,iRpc,iBin,dCorP,dCorP+dYShift);
1529 fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dYShift);
1530 fhCalChannelDy->Fill(dYShift);
1531 }
1532 }
1533 }
1534 else { // no individuak y-corrections
1535 for (Int_t iBin = 0; iBin < fhCorPos[iDetIndx]->GetNbinsX(); iBin++) {
1536 Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1537 fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dYShift);
1538 fhCalChannelDy->Fill(dYShift);
1539 }
1540 }
1541 } break; // case iOpt0=0
1542
1543 case 1: // update channel mean
1544 {
1545 LOG(info) << "Update time offsets for TSR " << iSmType << iSm << iRpc << " with opt1 " << iOpt1;
1546 TProfile* hpP = fhCalPos[iDetIndx]->ProfileX();
1547 TProfile* hpT = fhCalTOff[iDetIndx]->ProfileX();
1548 TH1* hCalP = fhCalPos[iDetIndx]->ProjectionX();
1549 TH1* hCalT = fhCalTOff[iDetIndx]->ProjectionX();
1550 //fhCorPos[iDetIndx]->Add((TH1 *)hpP,-1.);
1551 //fhCorTOff[iDetIndx]->Add((TH1*)hpT,-1.);
1552 if (iOpt1 == 3) { // update counter times from track pulls
1553 TH1* hCalTpy = fhCalTOff[iDetIndx]->ProjectionY();
1554 Double_t dCtsTpy = hCalTpy->GetEntries();
1555 Double_t dDt = 0;
1556 if (dCtsTpy > MINCTS) {
1557 double dFMean = hCalTpy->GetBinCenter(hCalTpy->GetMaximumBin());
1558 double dFLim = 0.5; // CAUTION, fixed numeric value
1559 double dBinSize = hCalTpy->GetBinWidth(1);
1560 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1561 TFitResultPtr fRes = hCalTpy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1562 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1563 dDt = fRes->Parameter(1);
1564 LOG(info) << "Update cor hist " << fhCorTOff[iDetIndx]->GetName() << " by " << dDt << ", " << dBeamTOff;
1565 if (iSmType == 5) { // do not shift beam counter in time
1566 dDt -= dBeamTOff;
1567 }
1568 else {
1569 dDt += dBeamTOff;
1570 }
1571 fhCalCounterDt->Fill(dDt);
1572 for (int iCh = 0; iCh < fhCorTOff[iDetIndx]->GetNbinsX(); iCh++) {
1573 fhCalChannelDt->Fill(dDt);
1574 fhCorTOff[iDetIndx]->SetBinContent(iCh + 1, fhCorTOff[iDetIndx]->GetBinContent(iCh + 1) + dDt);
1575 }
1576 }
1577 }
1578 }
1579 else {
1580 double dTOffMeanShift = 0.;
1581 double dNCh = 0.;
1582 for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
1583 Double_t dDt = hpT->GetBinContent(iBin + 1);
1584 Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
1585 Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
1586 Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
1587 Double_t dDp = hpP->GetBinContent(iBin + 1);
1588 Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1589 LOG(debug) << "Cts check for " << fhCalTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
1590 << ", " << dCtsP << ", " << MINCTS;
1591 if (dCtsT > MINCTS && dCtsP > MINCTS) {
1592 // Fit Gaussian around peak
1593 TH1* hpPy =
1594 (TH1*) fhCalPos[iDetIndx]->ProjectionY(Form("PosPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
1595 if (hpPy->Integral() > MINCTS) {
1596 Double_t dFMean = hpPy->GetBinCenter(hpPy->GetMaximumBin());
1597 Double_t dFLim = 0.5; // CAUTION, fixed numeric value
1598 Double_t dBinSize = hpPy->GetBinWidth(1);
1599 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1600 TFitResultPtr fRes = hpPy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1601 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1602 dDp = fRes->Parameter(1); //overwrite mean
1603 fhCalChannelDy->Fill(dDt);
1604 }
1605 }
1606 TH1* hpTy =
1607 (TH1*) fhCalTOff[iDetIndx]->ProjectionY(Form("TOffPy_%d_%d", iDetIndx, iBin), iBin + 1, iBin + 1);
1608 if (hpTy->Integral() > MINCTS) {
1609 Double_t dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
1610 Double_t dFLim = 0.5; // CAUTION, fixed numeric value
1611 Double_t dBinSize = hpTy->GetBinWidth(1);
1612 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1613 TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1614 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1615 dDt = fRes->Parameter(1); //overwrite mean
1616 fhCalChannelDt->Fill(dDt);
1617 }
1618 }
1619 // Double_t dDpRes = fRes->Parameter(2);
1620 if (iSmType == 5) { // do not shift beam counter in time
1621 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt - dBeamTOff);
1622 dTOffMeanShift += dDt - dBeamTOff;
1623 }
1624 else {
1625 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt + dBeamTOff);
1626 dTOffMeanShift += dDt + dBeamTOff;
1627 }
1628 dNCh += 1.;
1629 if (iOpt1 > 0) { // active y-position correction
1630 fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dDp);
1631 }
1632 if (iDetIndx > -1) {
1633 LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, beam %6.3f, new %6.3f",
1634 fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dBeamTOff,
1635 dCorT + dDt + dBeamTOff);
1636 }
1637 }
1638 }
1639 if (dNCh > 0) dTOffMeanShift /= dNCh;
1640 //LOG(info) << "Apply dTOffMeanShift " << dTOffMeanShift << " to " << fhCorTOff[iDetIndx]->GetName();
1641 for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) { //preserve mean offset
1642 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1,
1643 fhCorTOff[iDetIndx]->GetBinContent(iBin + 1) - dTOffMeanShift);
1644 }
1645 }
1646 } break;
1647 case 2: { // iopt1==11: update individual channel walks from inside Cluster deviations
1648 switch (iOpt1) {
1649 case 10: {
1650 LOG(debug) << "Update Cluster Offsets for TSR " << iSmType << iSm << iRpc;
1651 if (iSmType == 5) continue; // do not shift beam counter in time
1652 //TProfile* hpP = fhCalDelPos[iDetIndx]->ProfileX();
1653 TProfile* hpT = fhCalDelTOff[iDetIndx]->ProfileX();
1654 TH1* hCalP = fhCalDelPos[iDetIndx]->ProjectionX();
1655 TH1* hCalT = fhCalDelTOff[iDetIndx]->ProjectionX();
1656 //double dTOffMeanShift = 0.;
1657 //double dNCh = 0.;
1658 for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
1659 Double_t dDt = hpT->GetBinContent(iBin + 1);
1660 Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
1661 Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
1662 Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
1663 //Double_t dDp = hpP->GetBinContent(iBin + 1);
1664 //Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1665 LOG(debug) << "Cts check for " << fhCalDelTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
1666 << ", " << dCtsP << ", " << MINCTS;
1667 if (dCtsT > MINCTS) { //&& dCtsP > MINCTS) {
1668 // Fit Gaussian around peak
1669 TH1* hpTy = (TH1*) fhCalDelTOff[iDetIndx]->ProjectionY(Form("DelTOffPy_%d_%d", iDetIndx, iBin),
1670 iBin + 1, iBin + 1);
1671 double dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
1672 double dFLim = 0.5; // CAUTION, fixed numeric value
1673 double dBinSize = hpTy->GetBinWidth(1);
1674 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1675 TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1676 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1677 dDt = fRes->Parameter(1); //overwrite mean
1678 }
1679 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt);
1680 fhCalChannelDt->Fill(dDt);
1681 }
1682 if (iDetIndx > -1) {
1683 LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, new %6.3f",
1684 fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dCorT + dDt);
1685 }
1686 }
1687 } break;
1688 case 12: {
1689 LOG(debug) << "Update Cluster Position for TSR " << iSmType << iSm << iRpc;
1690 if (iSmType == 5) continue; // do not shift beam counter in time
1691 TProfile* hpP = fhCalDelPos[iDetIndx]->ProfileX();
1692 TProfile* hpT = fhCalDelTOff[iDetIndx]->ProfileX();
1693 TH1* hCalP = fhCalDelPos[iDetIndx]->ProjectionX();
1694 TH1* hCalT = fhCalDelTOff[iDetIndx]->ProjectionX();
1695 //double dTOffMeanShift = 0.;
1696 //double dNCh = 0.;
1697 for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
1698 Double_t dDt = hpT->GetBinContent(iBin + 1);
1699 Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
1700 Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
1701 Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
1702 Double_t dDp = hpP->GetBinContent(iBin + 1);
1703 Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1704 LOG(debug) << "Cts check for " << fhCalDelTOff[iDetIndx]->GetName() << ", bin " << iBin << ": " << dCtsT
1705 << ", " << dCtsP << ", " << MINCTS;
1706 if (dCtsT > MINCTS && dCtsP > MINCTS) {
1707 // Fit Gaussian around peak
1708 TH1* hpTy = (TH1*) fhCalDelTOff[iDetIndx]->ProjectionY(Form("DelTOffPy_%d_%d", iDetIndx, iBin),
1709 iBin + 1, iBin + 1);
1710 double dFMean = hpTy->GetBinCenter(hpTy->GetMaximumBin());
1711 double dFLim = 0.5; // CAUTION, fixed numeric value
1712 double dBinSize = hpTy->GetBinWidth(1);
1713 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1714 TFitResultPtr fRes = hpTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1715 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1716 dDt = fRes->Parameter(1); //overwrite mean
1717 }
1718 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dDt);
1719 fhCalChannelDt->Fill(dDt);
1720
1721 TH1* hpPy = (TH1*) fhCalDelPos[iDetIndx]->ProjectionY(Form("DelPosPy_%d_%d", iDetIndx, iBin),
1722 iBin + 1, iBin + 1);
1723 dFMean = hpPy->GetBinCenter(hpTy->GetMaximumBin());
1724 dFLim = 0.5; // CAUTION, fixed numeric value
1725 dBinSize = hpPy->GetBinWidth(1);
1726 dFLim = TMath::Max(dFLim, 5. * dBinSize);
1727 fRes = hpPy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1728 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1729 dDp = fRes->Parameter(1); //overwrite mean
1730 }
1731 fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dDp);
1732 fhCalChannelDy->Fill(dDp);
1733 }
1734 if (iDetIndx > -1) {
1735 LOG(debug) << Form("Update %s: bin %02d, Cts: %d, Old %6.3f, dev %6.3f, new %6.3f",
1736 fhCorTOff[iDetIndx]->GetName(), iBin, (Int_t) dCtsT, dCorT, dDt, dCorT + dDt);
1737 }
1738 }
1739 } break;
1740 case 0:
1741 case 11: {
1742 LOG(debug) << "Update Walks for TSR " << iSmType << iSm << iRpc;
1743 //if (iSmType == 5) continue; // no walk correction for beam counter up to now
1744 const Double_t MinCounts = 10.;
1745 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
1746 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
1747 TH1* hCY = fhCalPos[iDetIndx]->ProjectionY(Form("%s_py_%d", fhCalPos[iDetIndx]->GetName(), iCh),
1748 iCh + 1, iCh + 1);
1749 if (hCY->GetEntries() > 1) fhCalChannelDy->Fill(hCY->GetRMS());
1750 //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
1751 TProfile* hpW0 = fhCalWalk[iDetIndx][iCh][0]->ProfileX(); // mean deviation
1752 TH1* hCW0 = fhCalWalk[iDetIndx][iCh][0]->ProjectionX(); // contributing counts
1753 TProfile* hpW1 = fhCalWalk[iDetIndx][iCh][1]->ProfileX(); // mean deviation
1754 TH1* hCW1 = fhCalWalk[iDetIndx][iCh][1]->ProjectionX(); // contributing counts
1755 if (iOpt0 == 3) { // use average distributions
1756 hpW0 = fhCalWalkAv[iDetIndx]->ProfileX(); // mean deviation
1757 hCW0 = fhCalWalkAv[iDetIndx]->ProjectionX(); // contributing counts
1758 hpW1 = fhCalWalkAv[iDetIndx]->ProfileX(); // mean deviation
1759 hCW1 = fhCalWalkAv[iDetIndx]->ProjectionX(); // contributing counts
1760 }
1761 if (hCW0->GetEntries() == 0 || hCW1->GetEntries() == 0) {
1762 //LOG(info) << "No entries in " << hCW0->GetName();
1763 continue;
1764 }
1765 Double_t dCorT = 0;
1766 for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][0]->GetNbinsX(); iBin++) {
1767 Double_t dCts0 = hCW0->GetBinContent(iBin + 1);
1768 Double_t dCts1 = hCW1->GetBinContent(iBin + 1);
1769 Double_t dWOff0 = fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin + 1); // current value
1770 Double_t dWOff1 = fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin + 1); // current value
1771 if (iBin > 0 && dCts0 == 0 && dCts1 == 0) {
1772 fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1,
1773 fhCorWalk[iDetIndx][iCh][0]->GetBinContent(iBin));
1774 fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1,
1775 fhCorWalk[iDetIndx][iCh][1]->GetBinContent(iBin));
1776 }
1777 dCorT = 0.;
1778 if (dCts0 > MinCounts && dCts1 > MinCounts) {
1779 dCorT = 0.5 * (hpW0->GetBinContent(iBin + 1) + hpW1->GetBinContent(iBin + 1));
1780 fhCalChannelDt->Fill(dCorT);
1781 }
1782 fhCorWalk[iDetIndx][iCh][0]->SetBinContent(iBin + 1, dWOff0 + dCorT); //set new value
1783 fhCorWalk[iDetIndx][iCh][1]->SetBinContent(iBin + 1, dWOff1 + dCorT); //set new value
1784 if (iSmType == 0 && iSm == 0 && iRpc == 2) // debugging
1785 LOG(info) << "UpdWalk " << fhCorWalk[iDetIndx][iCh][0]->GetName() << " bin " << iBin << " Tot "
1786 << hCW0->GetBinCenter(iBin + 1) << ", " << hCW1->GetBinCenter(iBin + 1) << ": curVal "
1787 << dWOff0 << ", " << dWOff1 << ", Cts " << dCts0 << ", " << dCts1 << ", dCorT " << dCorT
1788 << ", newVal " << dWOff0 + dCorT << ", " << dWOff1 + dCorT;
1789 }
1790 // determine effective/count rate weighted mean
1791 /*
1792 Double_t dMean = 0;
1793 Double_t dCtsAll = 0.;
1794 for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
1795 Double_t dCts = hCW->GetBinContent(iBin + 1);
1796 Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1);
1797 if (dCts > MinCounts) {
1798 dCtsAll += dCts;
1799 dMean += dCts * dWOff;
1800 }
1801 }
1802 if (dCtsAll > 0.) dMean /= dCtsAll;
1803
1804 //LOG(info) << "Mean shift for TSRCS " << iSmType << iSm << iRpc << iCh << iSide << ": " << dMean;
1805 // keep mean value at 0
1806 for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++) {
1807 Double_t dWOff = fhCorWalk[iDetIndx][iCh][iSide]->GetBinContent(iBin + 1); // current value
1808 if (iSmType == 0 && iSm == 0 && iRpc == 2) // debugging
1809 LOG(info) << "UpdWalk " << fhCorWalk[iDetIndx][iCh][iSide]->GetName()
1810 << " bin " << iBin << " Tot " << hCW->GetBinCenter(iBin + 1)
1811 << ": curDev " << hpW->GetBinContent(iBin + 1)
1812 << ", Cts " << hCW->GetBinContent(iBin + 1)
1813 << ", all " << dCtsAll
1814 << ", oldWOff " << dWOff - hpW->GetBinContent(iBin + 1)
1815 << ", Mean " << dMean
1816 << ", newWOff " << dWOff - dMean;
1817 fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1, dWOff + dMean); //set new value
1818 }
1819 */
1820 }
1821 } break; //iopt1=11 end
1822 default:; LOG(info) << "Calibrator option " << iOpt << " not implemented ";
1823 } // switch iOpt1 end
1824 } break;
1825
1826 case 9: // update channel means on cluster level
1827 {
1828 // position, do Edge fit by default
1829 //LOG(info) << "Update Offsets for TSR " << iSmType << iSm << iRpc;
1830
1831 //TProfile* hpP = fhCalPosition[iDetIndx]->ProfileX();
1832 TH1* hCalP = fhCalPosition[iDetIndx]->ProjectionX();
1833
1834 for (Int_t iBin = 0; iBin < fhCorPos[iDetIndx]->GetNbinsX(); iBin++) {
1835 Double_t dCorP = fhCorPos[iDetIndx]->GetBinContent(iBin + 1);
1836 //Double_t dDp = hpP->GetBinContent(iBin + 1);
1837 Double_t dCtsP = hCalP->GetBinContent(iBin + 1);
1838 if (dCtsP > MINCTS) {
1839 TH1* hPos_py = fhCalPosition[iDetIndx]->ProjectionY(
1840 Form("%s_py_%d", fhCalPosition[iDetIndx]->GetName(), iBin), iBin + 1, iBin + 1);
1841 double dYShift = hPos_py->GetMean();
1842 int iChId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
1843 CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
1844 if (NULL == fChannelInfo) LOG(fatal) << Form("invalid ChannelInfo for 0x%08x", iChId);
1845
1846 double dThr = dCtsP / 100;
1847 if (dThr < 3.) dThr = 3.;
1848 double* dRes =
1849 fTofClusterizer->find_yedges((const char*) (hPos_py->GetName()), dThr, fChannelInfo->GetSizey());
1850 /*
1851 LOG(info) << Form("EdgeY for %s, TSR %d%d%d: DY %5.2f, Len %5.2f, Size %5.2f ",
1852 hPos_py->GetName(), iSmType, iSm, iRpc, dRes[1], dRes[0], fChannelInfo->GetSizey());
1853 */
1854 if (TMath::Abs(dRes[0] - fChannelInfo->GetSizey()) / fChannelInfo->GetSizey() < 0.1) {
1855 dYShift = dRes[1];
1856 }
1857 else {
1858 dYShift = hPos_py->GetMean();
1859 }
1860 // apply correction
1861 //LOG(info)<<Form("UpdateCalPos TSRC %d%d%d%2d: %6.2f -> %6.2f ",iSmType,iSm,iRpc,iBin,dCorP,dCorP+dYShift);
1862 fhCorPos[iDetIndx]->SetBinContent(iBin + 1, dCorP + dYShift);
1863 }
1864 }
1865
1866 //TofOff or TOff as input ?
1867 //TProfile* hpT = fhCalTofOff[iDetIndx]->ProfileX();
1868 //TH1* hCalT = fhCalTofOff[iDetIndx]->ProjectionX();
1869 TProfile* hpT = fhCalTOff[iDetIndx]->ProfileX();
1870 TH1* hCalT = fhCalTOff[iDetIndx]->ProjectionX();
1871 for (Int_t iBin = 0; iBin < fhCorTOff[iDetIndx]->GetNbinsX(); iBin++) {
1872 //Double_t dDt = hpT->GetBinContent(iBin + 1);
1873 Double_t dCorT = fhCorTOff[iDetIndx]->GetBinContent(iBin + 1);
1874 Double_t dCtsT = hCalT->GetBinContent(iBin + 1);
1875 if (dCtsT > MINCTS) {
1876 double dTmean = hpT->GetBinContent(iBin + 1);
1877 TH1* hTy = (TH1*) fhCalTofOff[iDetIndx]->ProjectionY(
1878 Form("%s_py%d", fhCalTofOff[iDetIndx]->GetName(), iBin), iBin + 1, iBin + 1);
1879 Double_t dNPeak = hTy->GetBinContent(hTy->GetMaximumBin());
1880 if (dNPeak > MINCTS * 0.1) {
1881 Double_t dFMean = hTy->GetBinCenter(hTy->GetMaximumBin());
1882 Double_t dFLim = 2.0; // CAUTION, fixed numeric value
1883 Double_t dBinSize = hTy->GetBinWidth(1);
1884 dFLim = TMath::Max(dFLim, 10. * dBinSize);
1885 TFitResultPtr fRes = hTy->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1886 //if (fRes == 0 )
1887 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1888 //if (TMath::Abs(TMean - fRes->Parameter(1)) > 1.)
1889 /*
1890 LOG(info) << "CalTOff "
1891 << Form("TSRC %d%d%d%d, %s gaus %8.2f %8.2f %8.2f for "
1892 "TM %8.2f, TBeam %6.2f",
1893 iSmType, iSm, iRpc, iBin, hTy->GetName(), fRes->Parameter(0), fRes->Parameter(1),
1894 fRes->Parameter(2), dTmean, dBeamTOff);
1895 */
1896 dTmean = fRes->Parameter(1); //overwrite mean
1897 }
1898 }
1899 /*
1900 LOG(info) << Form("UpdateCalTOff TSRC %d%d%d%2d, cts %d: %6.2f -> %6.2f, %6.2f ", iSmType, iSm, iRpc,
1901 iBin, (int) dCtsT, dCorT, dCorT + dTmean, dCorT + hpT->GetBinContent(iBin + 1));
1902 */
1903 fhCorTOff[iDetIndx]->SetBinContent(iBin + 1, dCorT + dTmean);
1904 }
1905 }
1906
1907 //Tot
1908 double dCalTotMean = fTofClusterizer->GetTotMean(); // Target value
1909 for (Int_t iBin = 0; iBin < fhCorTot[iDetIndx]->GetNbinsX(); iBin++) {
1910 int ib = iBin + 1;
1911 TH1* hbin = fhCalTot[iDetIndx]->ProjectionY(Form("bin%d", ib), ib, ib);
1912 double dTotMean = hbin->GetMean();
1913 // Do gaus fit around maximum
1914 Double_t dCtsTot = hbin->GetEntries();
1915 if (dCtsTot > MINCTS) {
1916 double dFMean = hbin->GetBinCenter(hbin->GetMaximumBin());
1917 double dFLim = dFMean * 0.5; // CAUTION,
1918 //LOG(info)<<Form("FitTot TSRC %d%d%d%2d: Mean %6.2f, Width %6.2f ",iSmType,iSm,iRpc,iBin,dFMean,dFLim);
1919 if (dFMean > 2.) {
1920 TFitResultPtr fRes = hbin->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
1921 if (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED")) {
1922 dTotMean = fRes->Parameter(1); //overwrite mean
1923 }
1924 else {
1925 LOG(warn) << "FitFail " << hbin->GetName();
1926 }
1927 }
1928 }
1929 else { //append to broken channel list
1930 int iCh = iBin % 2;
1931 int iSide = iBin - iCh * 2;
1932 CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, iCh, iSide);
1933 Int_t iChId = fTofId->SetDetectorInfo(xDetInfo);
1934 TString shcmd = Form("echo %d >> %s ", iChId, cBadChannelFile.Data());
1935 gSystem->Exec(shcmd.Data());
1936 }
1937 if (dTotMean > 0.) { // prevent zero-divide
1938 double dCorTot = fhCorTot[iDetIndx]->GetBinContent(ib);
1939 /*
1940 LOG(info)<<Form("UpdateCalTot TSRC %d%d%d%2d: %6.2f, %6.2f, %6.2f -> %6.2f ",iSmType,iSm,iRpc,iBin,
1941 dCorTot,dCalTotMean,dTotMean,dCorTot*dTotMean/dCalTotMean);
1942 */
1943 fhCorTot[iDetIndx]->SetBinContent(ib, dCorTot * dTotMean / dCalTotMean);
1944 }
1945 }
1946
1947 } break;
1948 default: LOG(fatal) << "No valid calibration mode " << iOpt0;
1949 } //switch( iOpt0) end
1950 }
1951 }
1952 else {
1953 // currently not needed
1954 }
1955 TString fFile = fCalParFile->GetName();
1956 if (!fFile.Contains("/")) {
1957 TFile* fCalParFileNew = new TFile(Form("New_%s", fCalParFile->GetName()), "RECREATE");
1958 WriteHist(fCalParFileNew);
1959 fCalParFileNew->Close();
1960 }
1961 else {
1962 WriteHist(fCalParFile);
1963 }
1964 fCalParFile->Close();
1965
1967 gFile = oldFile;
1968 //gDirectory = oldDir;
1969 gDirectory->cd(oldDir->GetPath());
1970
1971 return kTRUE;
1972}
1973
1974void CbmTofCalibrator::ReadHist(TFile* fHist)
1975{
1976 LOG(info) << "Read Cor histos from file " << fHist->GetName();
1977 if (0 == fhCorPos.size()) {
1978 Int_t iNbDet = fDigiBdfPar->GetNbDet();
1979 //LOG(info) << "resize histo vector for " << iNbDet << " detectors ";
1980 fhCorPos.resize(iNbDet);
1981 fhCorTOff.resize(iNbDet);
1982 fhCorTot.resize(iNbDet);
1983 fhCorTotOff.resize(iNbDet);
1984 fhCorSvel.resize(iNbDet);
1985 fhCorWalk.resize(iNbDet);
1986 }
1987
1988 for (Int_t iDetIndx = 0; iDetIndx < fDigiBdfPar->GetNbDet(); iDetIndx++) {
1989 Int_t iUniqueId = fDigiBdfPar->GetDetUId(iDetIndx);
1990 //Int_t iSmAddr = iUniqueId & DetMask;
1991 Int_t iSmType = CbmTofAddress::GetSmType(iUniqueId);
1992 Int_t iSm = CbmTofAddress::GetSmId(iUniqueId);
1993 Int_t iRpc = CbmTofAddress::GetRpcId(iUniqueId);
1994 //LOG(info) << "Get histo pointer for TSR " << iSmType<<iSm<<iRpc;
1995
1996 fhCorPos[iDetIndx] =
1997 (TH1*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
1998 if (NULL == fhCorPos[iDetIndx]) {
1999 LOG(error) << "hCorPos not found for TSR " << iSmType << iSm << iRpc;
2000 continue;
2001 }
2002 fhCorTOff[iDetIndx] =
2003 (TH1*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
2004 fhCorTot[iDetIndx] =
2005 (TH1*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
2006 fhCorTotOff[iDetIndx] =
2007 (TH1*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
2008
2009 Int_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
2010 fhCorWalk[iDetIndx].resize(iNbCh);
2011 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
2012 fhCorWalk[iDetIndx][iCh].resize(2);
2013 for (Int_t iSide = 0; iSide < 2; iSide++) {
2014 //LOG(info) << "Get walk histo pointer for TSRCS " << iSmType<<iSm<<iRpc<<iCh<<iSide;
2015 TString hname = Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_Walk_px", iSmType, iSm, iRpc, iCh, iSide);
2016 fhCorWalk[iDetIndx][iCh][iSide] = (TH1*) gDirectory->FindObjectAny(hname);
2017 if (NULL == fhCorWalk[iDetIndx][iCh][iSide]) {
2018 LOG(warn) << "No Walk histo for TSRCS " << iSmType << iSm << iRpc << iCh << iSide;
2019 if (NULL != fhCalWalk[iDetIndx][iCh][iSide]) {
2020 LOG(info) << "Create walk correction histo " << hname;
2021 fhCorWalk[iDetIndx][iCh][iSide] =
2022 new TH1D(hname,
2023 Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSm,
2024 iRpc, iCh, iSide),
2025 nbClWalkBinX, fTofClusterizer->GetTotMin(), fTofClusterizer->GetTotMax());
2026 }
2027 //continue;
2028 }
2029
2030 if (iSmType == 8 && iSide == 1) { //pad special treatment
2031 LOG(warn) << "Overwrite pad counter walk for TSRCS " << iSmType << iSm << iRpc << iCh << iSide;
2032 TH1* hTmp = (TH1*) gDirectory->FindObjectAny(
2033 Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%d_Walk_px", iSmType, iSm, iRpc, iCh, 1 - iSide));
2034 for (Int_t iBin = 0; iBin < fhCorWalk[iDetIndx][iCh][iSide]->GetNbinsX(); iBin++)
2035 fhCorWalk[iDetIndx][iCh][iSide]->SetBinContent(iBin + 1, hTmp->GetBinContent(iBin + 1));
2036 }
2037 }
2038 }
2039 }
2040
2041 int iNbTypes = fDigiBdfPar->GetNbSmTypes();
2042 hSvel.resize(iNbTypes);
2043 for (int iSmType = 0; iSmType < iNbTypes; iSmType++) {
2044 int iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
2045 if (iNbRpc > 0 && fDigiBdfPar->GetNbSm(iSmType) > 0) {
2046 TH1* hTmp = (TH1*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
2047 if (hTmp == NULL) {
2048 TDirectory* oldir = gDirectory;
2049 gROOT->cd();
2050 LOG(warn) << Form("cl_SmT%01d_Svel not found, creating ... in %s ", iSmType, gDirectory->GetName());
2051 hSvel[iSmType] =
2052 new TH1F(Form("cl_SmT%01d_Svel", iSmType), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iSmType),
2053 fDigiBdfPar->GetNbSm(iSmType) * iNbRpc, 0, fDigiBdfPar->GetNbSm(iSmType) * iNbRpc);
2054
2055 for (int i = 0; i < hSvel[iSmType]->GetNbinsX(); i++)
2056 hSvel[iSmType]->SetBinContent(i + 1, 1.); // set default
2057 oldir->cd();
2058 }
2059 else {
2060 //LOG(info) << "hSvel "<< hTmp->GetName() << " located in " << gDirectory->GetName();
2061 TDirectory* oldir = gDirectory;
2062 gROOT->cd();
2063 hSvel[iSmType] = (TProfile*) hTmp->Clone();
2064 //LOG(info) << "hSvel "<< hSvel[iSmType]->GetName() << " cloned in " << gDirectory->GetName();
2065 oldir->cd();
2066 }
2067 }
2068 }
2069}
2070
2071void CbmTofCalibrator::WriteHist(TFile* fHist)
2072{
2073 LOG(info) << "Write Cor histos to file " << fHist->GetName();
2074 TDirectory* oldir = gDirectory;
2075 fHist->cd();
2076 for (Int_t iDetIndx = 0; iDetIndx < fDigiBdfPar->GetNbDet(); iDetIndx++) {
2077 if (NULL == fhCorPos[iDetIndx]) continue;
2078 fhCorPos[iDetIndx]->Write();
2079 fhCorTOff[iDetIndx]->Write();
2080 fhCorTot[iDetIndx]->Write();
2081 fhCorTotOff[iDetIndx]->Write();
2082
2083 Int_t iNbCh = (Int_t) fhCorWalk[iDetIndx].size();
2084 for (Int_t iCh = 0; iCh < iNbCh; iCh++) {
2085 for (Int_t iSide = 0; iSide < 2; iSide++) {
2086 //LOG(info)<<"Write " << fhCorWalk[iDetIndx][iCh][iSide]->GetName();
2087 if (NULL != fhCorWalk[iDetIndx][iCh][iSide])
2088 fhCorWalk[iDetIndx][iCh][iSide]->Write();
2089 else
2090 LOG(warn) << "Walk correction histo for DetIndx " << iDetIndx << ", Ch " << iCh << ", S " << iSide
2091 << " not found";
2092 }
2093 }
2094 }
2095 for (Int_t iS = 0; iS < fDigiBdfPar->GetNbSmTypes(); iS++) {
2096 if (NULL != hSvel[iS]) {
2097 hSvel[iS]->Write();
2098 LOG(warn) << "Wrote " << hSvel[iS]->GetName() << " to " << gDirectory->GetName();
2099 }
2100 else {
2101 LOG(warn) << Form("cl_SmT%01d_Svel not found ", iS);
2102 }
2103 }
2104 oldir->cd();
2105}
2106
2107void CbmTofCalibrator::HstDoublets(CbmTofTracklet* pTrk)
2108{
2109 for (Int_t iHit = 0; iHit < pTrk->GetNofHits() - 1; iHit++) {
2110 for (Int_t iHit1 = 0; iHit1 < pTrk->GetNofHits(); iHit1++) {
2111 if (iHit == iHit1) continue;
2112 CbmTofHit* pHit0 = pTrk->GetTofHitPointer(iHit);
2113 CbmTofHit* pHit1 = pTrk->GetTofHitPointer(iHit1);
2114 //auto iHind0=pTrk->GetTofHitIndex(iHit);
2115 //auto iHind1=pTrk->GetTofHitIndex(iHit+1);
2116 int iDind0 = fDetIdIndexMap[pTrk->GetTofDetIndex(iHit)];
2117 int iDind1 = fDetIdIndexMap[pTrk->GetTofDetIndex(iHit1)];
2118 if (pHit1->GetZ() < pHit0->GetZ()) {
2119 continue;
2120 //iHind0=pTrk->GetTofHitIndex(iHit+1);
2121 //iHind1=pTrk->GetTofHitIndex(iHit);
2122 /*
2123 iDind0=fDetIdIndexMap[ pTrk->GetTofDetIndex(iHit1) ]; // numbering according to BDF
2124 iDind1=fDetIdIndexMap[ pTrk->GetTofDetIndex(iHit) ];
2125 pHit0=pTrk->GetTofHitPointer(iHit1);
2126 pHit1=pTrk->GetTofHitPointer(iHit);
2127 */
2128 }
2129 int iHst = iDind0 * 100 + iDind1;
2130 if (NULL == fhDoubletDt[iHst]) {
2131 // create histograms
2132 TDirectory* oldir =
2133 gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
2134 gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
2135 TString hNameDt = Form("hDoubletDt_%02d%02d", iDind0, iDind1);
2136 LOG(info) << Form("Book histo %lu %s, Addrs 0x%08x, 0x%08x ", fhDoubletDt.size(), hNameDt.Data(),
2137 pTrk->GetTofDetIndex(iHit), pTrk->GetTofDetIndex(iHit + 1));
2138 TH1D* pHstDt = new TH1D(hNameDt, Form("%s; #Delta t (ns); cts", hNameDt.Data()), 100, -5., 5.);
2139 fhDoubletDt[iHst] = pHstDt;
2140 TString hNameDd = Form("hDoubletDd_%02d%02d", iDind0, iDind1);
2141 TH1D* pHstDd = new TH1D(hNameDd, Form("%s; #Delta D (cm); cts", hNameDd.Data()), 200, 0., 200.);
2142 fhDoubletDd[iHst] = pHstDd;
2143 TString hNameV = Form("hDoubletV_%02d%02d", iDind0, iDind1);
2144 TH1D* pHstV = new TH1D(hNameV, Form("%s; v (cm/ns); cts", hNameV.Data()), 100, 0., 100.);
2145 fhDoubletV[iHst] = pHstV;
2146 oldir->cd();
2147 }
2148 // Fill Histograms
2149 double dDt = pHit1->GetTime() - pHit0->GetTime();
2150 double dDd = pTrk->Dist3D(pHit1, pHit0);
2151 fhDoubletDt[iHst]->Fill(dDt);
2152 fhDoubletDd[iHst]->Fill(dDd);
2153 fhDoubletV[iHst]->Fill(dDd / dDt);
2154 }
2155 }
2156}
2157
2158double* CbmTofCalibrator::find_tofedge(const char* hname, double dThr, double dLen)
2159{
2160 TH1* h1 = (TH1*) gROOT->FindObjectAny(hname);
2161 static double dRes[2] = {2 * 0};
2162 dRes[0] = 0.;
2163 if (NULL != h1) {
2164 /*
2165 LOG(debug) << "Inspect " << h1->GetName() << ", entries " << h1->GetNbinsX()
2166 << ", Thr " << dThr << ", Len " << dLen << endl;
2167 */
2168
2169 const int iMaxInt = 3;
2170 int iMax = h1->GetMaximumBin();
2171 int iLow = TMath::Max(iMax - iMaxInt, 1);
2172 int iHigh = TMath::Min(iMax + iMaxInt, h1->GetNbinsX());
2173 double dMax = 0.;
2174 int nMax = 0;
2175 for (int i = iLow; i < iHigh; i++) {
2176 dMax += h1->GetBinContent(i);
2177 nMax++;
2178 }
2179 dMax /= nMax;
2180
2181 // determine average noise and RMS per bin
2182 double dNoise = 0;
2183 double dNoise2 = 0;
2184 double dRMS = 0;
2185 int iNbins = 0;
2186 for (int iBin = 0; iBin < h1->GetNbinsX() / 4; iBin++) {
2187 if (h1->GetBinContent(iBin + 1) > 0 || iNbins > 0) { // start counting from first bin with entries
2188 dNoise += h1->GetBinContent(iBin + 1);
2189 dNoise2 += h1->GetBinContent(iBin + 1) * h1->GetBinContent(iBin + 1);
2190 iNbins++;
2191 }
2192 }
2193 if (iNbins > 0) {
2194 dNoise /= iNbins;
2195 dNoise2 /= iNbins;
2196 dRMS = TMath::Sqrt(dNoise2 - dNoise * dNoise);
2197 }
2198 else {
2199 dNoise = h1->GetMaximum() * 0.1;
2200 LOG(warn) << h1->GetName() << ": Noise level could not be determined, init to " << dNoise;
2201 }
2202 double dLev = dNoise + (dMax - dNoise) * dThr;
2203 dLev = TMath::Max(dNoise + 3. * dRMS, dLev);
2204 if (dLev == 0) {
2205 LOG(warn) << "Invalid threshold level " << dLev;
2206 return dRes;
2207 }
2208 int iBl = -1;
2209 double xLow = h1->GetBinCenter(1);
2210
2211 for (int iBin = 1; iBin < h1->GetNbinsX() - dLen; iBin++) {
2212 int iLen = 0;
2213 for (; iLen < dLen; iLen++) {
2214 LOG(debug) << h1->GetName() << ": Lev " << dLev << ", Bin " << iBin << ", Len " << iLen << ", cont "
2215 << h1->GetBinContent(iBin + iLen);
2216 if (h1->GetBinContent(iBin + iLen) < dLev) break;
2217 }
2218 if (iLen == dLen) {
2219 if (iBin == 1) {
2220 LOG(warn) << "find_tofedge: " << h1->GetName() << ", Lvl1 " << dLev;
2221 iBin--;
2222 dLev *= 2;
2223 continue;
2224 }
2225 iBl = iBin;
2226 xLow = h1->GetBinCenter(iBin);
2227 if (iBin == 1) { // not active since did not work, case already caught above
2228 LOG(warn) << GetName() << ": " << h1->GetName() << " Lvl1 " << dLev;
2229 xLow *= 2;
2230 }
2231 break;
2232 }
2233 }
2234 dRes[0] = xLow;
2235 LOG(info) << h1->GetName() << " with Thr " << dLev << ", Noise " << dNoise << ", RMS " << dRMS << " at iBl " << iBl
2236 << ", TOff " << dRes[0];
2237 }
2238 return dRes;
2239}
2240
2241double* CbmTofCalibrator::find_tofedge(const char* hname)
2242{
2243 double dThr = 0.1;
2244 double dLen = 3.;
2245 return find_tofedge(hname, dThr, dLen);
2246}
2247
2248double CbmTofCalibrator::CalcChi2(TH1* h1, TH1* h2, int i)
2249{
2250 double chi2 = 0.;
2251 double nBin = 0.;
2252 double nEntries1 = h1->GetEntries();
2253 double nEntries2 = h2->GetEntries();
2254 const int CompRange = (int) h1->GetNbinsX() / 4;
2255 for (int iBin1 = CompRange; iBin1 < h1->GetNbinsX() - CompRange; iBin1++) {
2256 double diff = h1->GetBinContent(iBin1) / nEntries1 - h2->GetBinContent(iBin1 + i) / nEntries2;
2257 chi2 += diff * diff;
2258 nBin++;
2259 }
2260 chi2 /= nBin;
2261 return chi2;
2262}
2263
2264
2265Double_t CbmTofCalibrator::f1_tedge(double* x, double* par)
2266{
2267 double xx = x[0];
2268 //double wx = par[0] + par[1] * TMath::Exp(xx*par[2]);
2269 double wx = par[0] + par[1] * (1. + TMath::Erf((xx - par[2]) / par[3]));
2270
2271 return wx;
2272}
2273
2274double* CbmTofCalibrator::fit_tofedge(const char* hname, Double_t TMax, Double_t dThr)
2275{
2276 TH1* h1;
2277 static double res[10] = {10 * 0.};
2278 double err[10];
2279 h1 = (TH1*) gROOT->FindObjectAny(hname);
2280 if (NULL != h1) {
2281 const double MinCts = 1000.;
2282 res[2] = 0.;
2283 if (h1->GetEntries() < MinCts) {
2284 LOG(warn) << h1->GetName() << ": too few entries for edgefit " << h1->GetEntries();
2285 return &res[2];
2286 }
2287
2288 TAxis* xaxis = h1->GetXaxis();
2289 Double_t Tmin = xaxis->GetXmin();
2290 Double_t Tmax = xaxis->GetXmax();
2291 Double_t dXmax = Tmax;
2292 if (TMax == 0.)
2293 TMax = Tmax;
2294 else
2295 Tmax = TMax - 1.;
2296 Tmax = TMath::Max(Tmax, Tmin + 0.3);
2297 TFitResultPtr fRes = h1->Fit("pol0", "SQM", "", Tmin, Tmax);
2298 TF1* f0 = h1->GetFunction("pol0");
2299 if (fRes != 0 && f0 != NULL && (gMinuit->fCstatu.Contains("OK") || gMinuit->fCstatu.Contains("CONVERGED"))) {
2300 //double dBackground=fRes->Parameter(0) ;
2301 double dBackground = f0->GetParameter(0);
2302 //cout << h1->GetName() << ": av background rate \t" << dBackground<< endl;
2303 TF1* f1 = new TF1("TEdge", f1_tedge, Tmin, dXmax, 4);
2304 if (TMax != 0.) {
2305 f1->SetParameters(dBackground, 1000., 0., 0.1);
2306 f1->SetParLimits(0, dBackground, dBackground); //fix background
2307 }
2308 f1->SetParLimits(2, -0.5, 0.5);
2309 f1->SetParLimits(3, 0.05, 0.5);
2310
2311 h1->Fit("TEdge", "SQM", "", Tmin, TMax);
2312 res[9] = f1->GetChisquare();
2313
2314 for (int i = 0; i < 4; i++) {
2315 res[i] = f1->GetParameter(i);
2316 err[i] = f1->GetParError(i);
2317 }
2318
2319 if (res[1] < 0.) {
2320 LOG(warn) << h1->GetName() << "fitted wrong sign of step function ";
2321 res[0] = h1->GetBinCenter(h1->GetNbinsX());
2322 return &res[0];
2323 }
2324 // determine toff at 10*background
2325 /*
2326 double dToff=0.;
2327 if (res[1] !=0. && res[2] != 0.) dToff=TMath::Log(10*dBackground/res[1])/res[2];
2328 LOG(info) << "TEdge Fit of " << hname
2329 << " in [" << Tmin <<","<< TMax <<"]"
2330 << " ended with chi2 = " << res[9]
2331 << ", Toff " << dToff
2332 << Form(", noise level %7.2f +/- %5.2f, e0 %7.2f "
2333 "+/- %5.2f, exp = %7.2f +/- %5.2f",
2334 res[0], err[0], res[1], err[1], res[2], err[2]);
2335 res[0]=dToff;
2336 */
2337 // FIXME: code duplication from find_tofedge
2338 const int iMaxInt = 3;
2339 int iMax = h1->GetMaximumBin();
2340 int iLow = TMath::Max(iMax - iMaxInt, 1);
2341 int iHigh = TMath::Min(iMax + iMaxInt, h1->GetNbinsX());
2342 double dMax = 0.;
2343 int nMax = 0;
2344 for (int i = iLow; i < iHigh; i++) {
2345 dMax += h1->GetBinContent(i);
2346 nMax++;
2347 }
2348 dMax /= nMax;
2349 double dLev = dThr * dMax;
2350 dLev = dBackground + (dMax - dBackground) * 0.1;
2351
2352 double dToff = f1->GetX(dLev, Tmin, dXmax);
2353
2354 LOG(info) << "TEdge Fit of " << hname << " in [" << Tmin << "," << TMax << "/" << dXmax << "]"
2355 << ": chi2 " << res[9]
2356 << Form(
2357 ", noise %6.2f +/-%4.2f, norm %6.2f +/-%4.2f, mean %6.2f +/-%5.3f, wid %6.2f +/-%5.3f -> %6.3f ",
2358 res[0], err[0], res[1], err[1], res[2], err[2], res[3], err[3], dToff);
2359 res[2] = dToff;
2360 }
2361 else {
2362 LOG(warn) << "pol0 fit failed, Minuit " << gMinuit->fCstatu << ", fRes " << fRes << ", "
2363 << ", f0 " << f0;
2364 }
2365 }
2366 return &res[2];
2367}
2368
2369double* CbmTofCalibrator::fit_tofedge(const char* hname)
2370{
2371 Double_t Tmax = 0.3;
2372 return fit_tofedge(hname, Tmax, fTofClusterizer->GetEdgeThr());
2373}
2374
2375double CbmTofCalibrator::TruncatedMeanY(TH2* h2, double dRmsLim)
2376{
2377 double dMeanIn = h2->ProjectionY()->GetMean();
2378 double dRmsIn = h2->ProjectionY()->GetRMS();
2379 const int Nx = h2->GetNbinsX();
2380
2381 int Ntru = 0;
2382 TH1* h1[Nx];
2383 double dMean = 0.;
2384 double dMean2 = 0.;
2385 double dNorm = 0.;
2386 double dRms = 0.;
2387 for (int i = 0; i < Nx; i++) {
2388 h1[i] = h2->ProjectionY(Form("%s_py%d", h2->GetName(), i), i + 1, i + 1);
2389 if (TMath::Abs(dMeanIn - h1[i]->GetMean()) < dRmsLim * dRmsIn) {
2390 dMean += h1[i]->GetMean() * h1[i]->GetEntries();
2391 dMean2 += h1[i]->GetMean() * h1[i]->GetMean() * h1[i]->GetEntries();
2392 dNorm += h1[i]->GetEntries();
2393 Ntru++;
2394 }
2395 }
2396 if (Ntru > 0) {
2397 dMean /= dNorm;
2398 dMean2 /= dNorm;
2399 if (dMean2 > dMean * dMean) {
2400 dRms = TMath::Sqrt(dMean2 - dMean * dMean);
2401 }
2402 else {
2403 LOG(error) << h2->GetName() << " Invalid Rms calculation: " << dMean2 << ", " << dMean * dMean;
2404 return dMean;
2405 }
2406 double dMean1 = dMean;
2407 dMean = 0.;
2408 dMean2 = 0.;
2409 dNorm = 0.;
2410 Ntru = 0;
2411 for (int j = 0; j < Nx; j++) {
2412 if (TMath::Abs(dMean1 - h1[j]->GetMean()) < dRmsLim * dRms) {
2413 dMean += h1[j]->GetMean() * h1[j]->GetEntries();
2414 dMean2 += h1[j]->GetMean() * h1[j]->GetMean() * h1[j]->GetEntries();
2415 dNorm += h1[j]->GetEntries();
2416 Ntru++;
2417 }
2418 }
2419 if (Ntru > 0) {
2420 dMean /= dNorm;
2421 }
2422 }
2423 LOG(info) << h2->GetName() << ": TruncY " << dMean << "( " << Ntru << ")"
2424 << ", full " << dMeanIn << ", RMS " << dRms;
2425 return dMean;
2426}
2427
2428ClassImp(CbmTofCalibrator)
ClassImp
ClassImp(CbmConverterManager)
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
pRef
CbmTofDigi * pRef
Definition CbmDeviceHitBuilderTof.cxx:86
dTmax
static Double_t dTmax
Definition CbmDeviceHitBuilderTof.cxx:84
fHist
static TFile * fHist
Definition CbmHadronAnalysis.cxx:82
DetMask
const Int_t DetMask
Definition CbmTofAnaTestbeam.cxx:79
fDigiMan
CbmDigiManager * fDigiMan
Definition CbmTofAnaTestbeam.cxx:91
hSvel
std::vector< TH1 * > hSvel
Definition CbmTofCalibrator.cxx:52
NevtT
static Int_t NevtT
Definition CbmTofCalibrator.cxx:760
dValidEdge
const double dValidEdge
Definition CbmTofCalibrator.cxx:50
MaxShift
const int MaxShift
Definition CbmTofCalibrator.cxx:47
cBadChannelFile
const TString cBadChannelFile
Definition CbmTofCalibrator.cxx:49
NevtH
static Int_t NevtH
Definition CbmTofCalibrator.cxx:358
cLight
const double cLight
Definition CbmTofCalibrator.cxx:48
nbClWalkBinX
const Int_t nbClWalkBinX
Definition CbmTofClusterizersDef.h:52
nbClWalkBinY
const Int_t nbClWalkBinY
Definition CbmTofClusterizersDef.h:53
CbmTofDetectorId_v21a.h
size
static constexpr size_t size()
Definition KfSimdPseudo.h:2
CbmDigiManager::IsPresent
static Bool_t IsPresent(ECbmModuleId systemId)
Presence of a digi branch.
Definition CbmDigiManager.cxx:121
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition CbmDigiManager.cxx:77
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition CbmDigiManager.h:126
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
CbmEvent::GetNofData
size_t GetNofData() const
Definition CbmEvent.cxx:53
CbmEvent::GetIndex
uint32_t GetIndex(ECbmDataType type, uint32_t iData)
Definition CbmEvent.cxx:42
CbmHit::GetTime
double GetTime() const
Definition CbmHit.h:76
CbmHit::GetAddress
int32_t GetAddress() const
Definition CbmHit.h:74
CbmHit::GetZ
double GetZ() const
Definition CbmHit.h:71
CbmMatch::GetLink
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmTofAddress::GetUniqueAddress
static uint32_t GetUniqueAddress(uint32_t Sm, uint32_t Rpc, uint32_t Channel, uint32_t Side=0, uint32_t SmType=0, uint32_t RpcType=0)
Definition CbmTofAddress.h:135
CbmTofAddress::GetSmId
static int32_t GetSmId(uint32_t address)
Definition CbmTofAddress.h:83
CbmTofAddress::GetRpcId
static int32_t GetRpcId(uint32_t address)
Definition CbmTofAddress.h:93
CbmTofAddress::GetSmType
static int32_t GetSmType(uint32_t address)
Definition CbmTofAddress.h:88
CbmTofAddress::GetChannelId
static int32_t GetChannelId(uint32_t address)
Definition CbmTofAddress.h:98
CbmTofCalibrator
contains filling and updating of calibration histos
Definition CbmTofCalibrator.h:51
CbmTofCalibrator::fhCalXYTOff
std::vector< TH3 * > fhCalXYTOff
Definition CbmTofCalibrator.h:118
CbmTofCalibrator::fhCalCluSize
std::vector< TH2 * > fhCalCluSize
Definition CbmTofCalibrator.h:114
CbmTofCalibrator::FillCalHist
void FillCalHist(CbmTofHit *pHit, Int_t iOpt, CbmEvent *pEvent=NULL)
Definition CbmTofCalibrator.cxx:359
CbmTofCalibrator::fhCalXYTot
std::vector< TH3 * > fhCalXYTot
Definition CbmTofCalibrator.h:119
CbmTofCalibrator::fTofCalDigiVec
const std::vector< CbmTofDigi > * fTofCalDigiVec
Definition CbmTofCalibrator.h:87
CbmTofCalibrator::fTofFindTracks
CbmTofFindTracks * fTofFindTracks
Definition CbmTofCalibrator.h:84
CbmTofCalibrator::fhCorTot
std::vector< TH1 * > fhCorTot
Definition CbmTofCalibrator.h:125
CbmTofCalibrator::fhDoubletDt
std::map< int, TH1 * > fhDoubletDt
Definition CbmTofCalibrator.h:131
CbmTofCalibrator::fTrackletTools
CbmTofTrackletTools * fTrackletTools
Definition CbmTofCalibrator.h:85
CbmTofCalibrator::fhCalPosition
std::vector< TH2 * > fhCalPosition
Definition CbmTofCalibrator.h:107
CbmTofCalibrator::fhCalWalkAv
std::vector< TH2 * > fhCalWalkAv
Definition CbmTofCalibrator.h:115
CbmTofCalibrator::fhCalDelTOff
std::vector< TH2 * > fhCalDelTOff
Definition CbmTofCalibrator.h:112
CbmTofCalibrator::fhCalTOff
std::vector< TH2 * > fhCalTOff
Definition CbmTofCalibrator.h:109
CbmTofCalibrator::f1_tedge
static double f1_tedge(double *x, double *par)
Definition CbmTofCalibrator.cxx:2265
CbmTofCalibrator::fhCalTotYWalk
std::vector< std::vector< std::vector< TH3 * > > > fhCalTotYWalk
Definition CbmTofCalibrator.h:120
CbmTofCalibrator::fhCalPos
std::vector< TH2 * > fhCalPos
Definition CbmTofCalibrator.h:108
CbmTofCalibrator::fEvtHeader
FairEventHeader * fEvtHeader
Definition CbmTofCalibrator.h:95
CbmTofCalibrator::fhCorSvel
std::vector< TH1 * > fhCorSvel
Definition CbmTofCalibrator.h:127
CbmTofCalibrator::fhCalDtWalk
std::vector< std::vector< std::vector< TH2 * > > > fhCalDtWalk
Definition CbmTofCalibrator.h:117
CbmTofCalibrator::ReadHist
void ReadHist(TFile *fhFile)
Definition CbmTofCalibrator.cxx:1974
CbmTofCalibrator::fTofId
CbmTofDetectorId * fTofId
Definition CbmTofCalibrator.h:92
CbmTofCalibrator::fhCalTofOff
std::vector< TH2 * > fhCalTofOff
Definition CbmTofCalibrator.h:110
CbmTofCalibrator::fhDoubletV
std::map< int, TH1 * > fhDoubletV
Definition CbmTofCalibrator.h:133
CbmTofCalibrator::fDigiBdfPar
CbmTofDigiBdfPar * fDigiBdfPar
Definition CbmTofCalibrator.h:90
CbmTofCalibrator::fhCalTot
std::vector< TH2 * > fhCalTot
Definition CbmTofCalibrator.h:106
CbmTofCalibrator::find_tofedge
double * find_tofedge(const char *hname, Double_t dThr, Double_t dLen)
CbmTofCalibrator::fTofDigiMatchColl
TClonesArray * fTofDigiMatchColl
Definition CbmTofCalibrator.h:93
CbmTofCalibrator::fhCorPos
std::vector< TH1 * > fhCorPos
Definition CbmTofCalibrator.h:123
CbmTofCalibrator::fit_tofedge
double * fit_tofedge(const char *hname, Double_t dTmax, Double_t dThr)
Definition CbmTofCalibrator.cxx:2274
CbmTofCalibrator::WriteHist
void WriteHist(TFile *fhFile)
Definition CbmTofCalibrator.cxx:2071
CbmTofCalibrator::fDetIdIndexMap
std::map< UInt_t, UInt_t > fDetIdIndexMap
Definition CbmTofCalibrator.h:130
CbmTofCalibrator::FillHitCalHist
void FillHitCalHist(CbmTofHit *pHit, Int_t iOpt, CbmEvent *pEvent=NULL, TClonesArray *tHitColl=NULL)
Definition CbmTofCalibrator.cxx:504
CbmTofCalibrator::HstDoublets
void HstDoublets(CbmTofTracklet *pTrk)
Definition CbmTofCalibrator.cxx:2107
CbmTofCalibrator::UpdateCalHist
Bool_t UpdateCalHist(Int_t iOpt)
Definition CbmTofCalibrator.cxx:982
CbmTofCalibrator::CalcChi2
double CalcChi2(TH1 *h1, TH1 *h2, int iShift)
Definition CbmTofCalibrator.cxx:2248
CbmTofCalibrator::fhCorTOff
std::vector< TH1 * > fhCorTOff
Definition CbmTofCalibrator.h:124
CbmTofCalibrator::fhCalCluTrms
std::vector< TH2 * > fhCalCluTrms
Definition CbmTofCalibrator.h:113
CbmTofCalibrator::TruncatedMeanY
double TruncatedMeanY(TH2 *pHst, double RmsLim=1.)
Definition CbmTofCalibrator.cxx:2375
CbmTofCalibrator::fhCorWalk
std::vector< std::vector< std::vector< TH1 * > > > fhCorWalk
Definition CbmTofCalibrator.h:128
CbmTofCalibrator::fhCorTotOff
std::vector< TH1 * > fhCorTotOff
Definition CbmTofCalibrator.h:126
CbmTofCalibrator::fhCalTotYTOff
std::vector< std::vector< std::vector< TH3 * > > > fhCalTotYTOff
Definition CbmTofCalibrator.h:121
CbmTofCalibrator::fhCalDelPos
std::vector< TH2 * > fhCalDelPos
Definition CbmTofCalibrator.h:111
CbmTofCalibrator::fTofClusterizer
CbmTofEventClusterizer * fTofClusterizer
Definition CbmTofCalibrator.h:83
CbmTofCalibrator::fhCalWalk
std::vector< std::vector< std::vector< TH2 * > > > fhCalWalk
Definition CbmTofCalibrator.h:116
CbmTofCalibrator::fDigiPar
CbmTofDigiPar * fDigiPar
Definition CbmTofCalibrator.h:89
CbmTofCalibrator::fhDoubletDd
std::map< int, TH1 * > fhDoubletDd
Definition CbmTofCalibrator.h:132
CbmTofCalibrator::fDigiMan
CbmDigiManager * fDigiMan
Definition CbmTofCalibrator.h:82
CbmTofCell::GetSizey
Double_t GetSizey() const
Definition CbmTofCell.h:40
CbmTofDetectorId_v21a
Definition CbmTofDetectorId_v21a.h:41
CbmTofDetectorId::SetDetectorInfo
virtual int32_t SetDetectorInfo(const CbmTofDetectorInfo detectorInfo)=0
CbmTofDigiBdfPar
Parameters class for the CBM ToF digitizer using beam data distributions.
Definition CbmTofDigiBdfPar.h:37
CbmTofDigiBdfPar::GetNbSmTypes
Int_t GetNbSmTypes() const
Definition CbmTofDigiBdfPar.h:62
CbmTofDigiBdfPar::GetNbSm
Int_t GetNbSm(Int_t iSmType) const
Definition CbmTofDigiBdfPar.cxx:528
CbmTofDigiBdfPar::GetNbChan
Int_t GetNbChan(Int_t iSmType, Int_t iRpc) const
Definition CbmTofDigiBdfPar.cxx:599
CbmTofDigiBdfPar::GetNbRpc
Int_t GetNbRpc(Int_t iSmType) const
Definition CbmTofDigiBdfPar.cxx:535
CbmTofDigiBdfPar::GetDetUId
Int_t GetDetUId(Int_t iDet)
Definition CbmTofDigiBdfPar.cxx:913
CbmTofDigiBdfPar::GetSigVel
Double_t GetSigVel(Int_t iSmType, Int_t iSm, Int_t iRpc) const
Definition CbmTofDigiBdfPar.cxx:564
CbmTofDigiPar::GetCell
CbmTofCell * GetCell(Int_t i)
Definition CbmTofDigiPar.h:50
CbmTofDigi
Data class for expanded digital TOF information.
Definition CbmTofDigi.h:47
CbmTofDigi::GetSide
double GetSide() const
Channel Side.
Definition CbmTofDigi.h:160
CbmTofDigi::GetChannel
double GetChannel() const
Channel .
Definition CbmTofDigi.h:156
CbmTofDigi::GetAddress
int32_t GetAddress() const
Inherited from CbmDigi.
Definition CbmTofDigi.h:112
CbmTofDigi::GetSm
double GetSm() const
Sm.
Definition CbmTofDigi.h:144
CbmTofDigi::GetTime
double GetTime() const
Inherited from CbmDigi.
Definition CbmTofDigi.h:131
CbmTofDigi::GetType
double GetType() const
Sm Type .
Definition CbmTofDigi.h:148
CbmTofDigi::GetRpc
double GetRpc() const
Detector aka Module aka RPC .
Definition CbmTofDigi.h:152
CbmTofDigi::GetTot
double GetTot() const
Alias for GetCharge.
Definition CbmTofDigi.h:140
CbmTofEventClusterizer::MasterToLocal
void MasterToLocal(const Int_t CellId, const Double_t *master, Double_t *local)
Definition CbmTofEventClusterizer.cxx:6757
CbmTofEventClusterizer::Instance
static CbmTofEventClusterizer * Instance()
Definition CbmTofEventClusterizer.h:69
CbmTofEventClusterizer::fit_ybox
virtual void fit_ybox(const char *hname)
Definition CbmTofEventClusterizer.cxx:5281
CbmTofEventClusterizer::GetHitPointer
CbmTofHit * GetHitPointer(int iHit)
Definition CbmTofEventClusterizer.cxx:6690
CbmTofEventClusterizer::GetMatchPointer
CbmMatch * GetMatchPointer(int iHit)
Definition CbmTofEventClusterizer.cxx:6698
CbmTofEventClusterizer::GetMatchIndexPointer
CbmMatch * GetMatchIndexPointer(int idx)
Definition CbmTofEventClusterizer.cxx:6707
CbmTofEventClusterizer::find_yedges
virtual double * find_yedges(const char *hname, Double_t dThr, Double_t dLen)
Definition CbmTofEventClusterizer.cxx:6635
CbmTofEventClusterizer::GetLocalY
double GetLocalY(CbmTofHit *pHit)
Definition CbmTofEventClusterizer.cxx:6736
CbmTofFindTracks::Instance
static CbmTofFindTracks * Instance()
Definition CbmTofFindTracks.h:69
CbmTofFindTracks::GetTofHitIndex
Int_t GetTofHitIndex(Int_t iHit)
Definition CbmTofFindTracks.h:187
CbmTofFindTracks::GetTtLight
Double_t GetTtLight() const
Definition CbmTofFindTracks.h:132
CbmTofFindTracks::GetTtTarg
Double_t GetTtTarg() const
Definition CbmTofFindTracks.h:131
CbmTofFindTracks::GetTOff
Double_t GetTOff(Int_t iAddr)
Definition CbmTofFindTracks.cxx:2237
CbmTofFindTracks::GetNReqStations
Int_t GetNReqStations() const
Definition CbmTofFindTracks.h:111
CbmTofHit::GetR
double GetR() const
Definition CbmTofHit.h:78
CbmTofHit::Dist3D
double Dist3D(CbmTofHit *pHit)
Definition CbmTofHit.h:85
CbmTofTrackletTools
contains fits and resolution functions
Definition CbmTofTrackletTools.h:25
CbmTofTrackletTools::GetTexpected
virtual Double_t GetTexpected(CbmTofTracklet *pTrk, Int_t iDetId, CbmTofHit *pHit, double TtLight=0.)
Definition CbmTofTrackletTools.cxx:265
CbmTofTracklet
Provides information on attaching a TofHit to a TofTrack.
Definition CbmTofTracklet.h:30
CbmTofTracklet::GetFitY
double GetFitY(double Z)
Definition CbmTofTracklet.cxx:505
CbmTofTracklet::GetTofHitIndex
int32_t GetTofHitIndex(int32_t ind) const
Definition CbmTofTracklet.h:72
CbmTofTracklet::ContainsAddr
virtual bool ContainsAddr(int32_t iAddr)
Definition CbmTofTracklet.cxx:450
CbmTofTracklet::GetTt
double GetTt() const
Definition CbmTofTracklet.h:59
CbmTofTracklet::GetNofHits
int32_t GetNofHits() const
Definition CbmTofTracklet.h:53
CbmTofTracklet::GetTofHitPointer
CbmTofHit * GetTofHitPointer(int32_t ind)
Definition CbmTofTracklet.h:74
CbmTofTracklet::GetFitT
double GetFitT(double R)
Definition CbmTofTracklet.cxx:507
CbmTofTracklet::GetTofDetIndex
int32_t GetTofDetIndex(int32_t ind) const
Definition CbmTofTracklet.h:75
CbmTofTracklet::GetFitX
double GetFitX(double Z)
Definition CbmTofTracklet.cxx:503
CbmTofTracklet::Dist3D
virtual double Dist3D(CbmTofHit *pHit0, CbmTofHit *pHit1)
Definition CbmTofTracklet.cxx:535