From bd90bf1da5e11800bb3733c8d0fc5948cffc36cd Mon Sep 17 00:00:00 2001
From: Norbert Herrmann <n.herrmann@physi.uni-heidelberg.de>
Date: Tue, 6 Apr 2021 15:32:41 +0200
Subject: [PATCH] add error estimate to Tof Hits

---
 macro/beamtime/pl_over_trk.C                  | 74 +++++++++++++------
 macro/beamtime/pl_pull_trk.C                  | 38 ++++------
 reco/detectors/tof/CbmTofCalibrator.cxx       |  4 +-
 reco/detectors/tof/CbmTofEventClusterizer.cxx |  9 ++-
 reco/detectors/tof/CbmTofFindTracks.cxx       | 22 +++---
 5 files changed, 86 insertions(+), 61 deletions(-)

diff --git a/macro/beamtime/pl_over_trk.C b/macro/beamtime/pl_over_trk.C
index 0524c9c935..b383955457 100644
--- a/macro/beamtime/pl_over_trk.C
+++ b/macro/beamtime/pl_over_trk.C
@@ -1,4 +1,5 @@
-void pl_over_trk(Int_t NSt = 4) {
+void pl_over_trk(Int_t NSt = 4)
+{
   //  TCanvas *can = new TCanvas("can22","can22");
   //  can->Divide(2,2);
   TCanvas* can = new TCanvas("can", "can", 50, 0, 800, 800);
@@ -26,7 +27,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -37,7 +39,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -48,7 +51,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -59,7 +63,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -74,12 +79,14 @@ void pl_over_trk(Int_t NSt = 4) {
         h1->Draw("");
         h1->SetMinimum(0.5);
         gPad->SetLogy();
-      } else {
+      }
+      else {
         h1->Draw("same");
       }
       if (iCol == 5 || iCol == 10) iCol++;
       h1->SetLineColor(iCol++);
-    } else {
+    }
+    else {
       cout << hname << " not found" << endl;
     }
   }
@@ -95,12 +102,14 @@ void pl_over_trk(Int_t NSt = 4) {
         h1->Draw("");
         h1->SetMinimum(0.5);
         gPad->SetLogy();
-      } else {
+      }
+      else {
         h1->Draw("same");
       }
       if (iCol == 5 || iCol == 10) iCol++;
       h1->SetLineColor(iCol++);
-    } else {
+    }
+    else {
       cout << hname << " not found" << endl;
     }
   }
@@ -116,12 +125,14 @@ void pl_over_trk(Int_t NSt = 4) {
         h1->Draw("");
         h1->SetMinimum(0.5);
         gPad->SetLogy();
-      } else {
+      }
+      else {
         h1->Draw("same");
       }
       if (iCol == 5 || iCol == 10) iCol++;
       h1->SetLineColor(iCol++);
-    } else {
+    }
+    else {
       cout << hname << " not found" << endl;
     }
   }
@@ -137,12 +148,14 @@ void pl_over_trk(Int_t NSt = 4) {
         h1->Draw("");
         h1->SetMinimum(0.5);
         gPad->SetLogy();
-      } else {
+      }
+      else {
         h1->Draw("same");
       }
       if (iCol == 5 || iCol == 10) iCol++;
       h1->SetLineColor(iCol++);
-    } else {
+    }
+    else {
       cout << hname << " not found" << endl;
     }
   }
@@ -154,7 +167,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -165,7 +179,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -176,7 +191,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -187,7 +203,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -206,13 +223,20 @@ void pl_over_trk(Int_t NSt = 4) {
   hname = Form("hCalDX0");
   h1    = (TH1*) gROOT->FindObjectAny(hname);
   if (h1 != NULL) {
-    h1->Draw("");
-    gPad->SetLogy();
     hname    = Form("hCalDY0");
     TH1* h1y = (TH1*) gROOT->FindObjectAny(hname);
+
+    Double_t dYmax = TMath::Max(h1->GetMaximum(), h1y->GetMaximum());
+    dYmax *= 1.1;
+    h1->SetMaximum(dYmax);
+
+    h1->Draw("");
+    gPad->SetLogy();
+
     h1y->SetLineColor(2);
     h1y->Draw("same");
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
 
@@ -225,7 +249,8 @@ void pl_over_trk(Int_t NSt = 4) {
     h1->SetMinimum(0.1 * h1->GetMaximum());
     gPad->SetLogy();
     hAll = (TH1*) h1->Clone();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
   hname = Form("hUsedHitsStation");
@@ -238,8 +263,8 @@ void pl_over_trk(Int_t NSt = 4) {
 
     can->cd(15);
     hEff->Draw();
-
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
   can->cd(16);
@@ -249,7 +274,8 @@ void pl_over_trk(Int_t NSt = 4) {
   if (h2 != NULL) {
     h2->Draw("colz");
     gPad->SetLogz();
-  } else {
+  }
+  else {
     cout << hname << " not found" << endl;
   }
   can->SaveAs("pl_over_trk.pdf");
diff --git a/macro/beamtime/pl_pull_trk.C b/macro/beamtime/pl_pull_trk.C
index d18f38f1b3..592a61da17 100644
--- a/macro/beamtime/pl_pull_trk.C
+++ b/macro/beamtime/pl_pull_trk.C
@@ -1,7 +1,5 @@
-void pl_pull_trk(Int_t NSt   = 8,
-                 Int_t iVar  = 0,
-                 Int_t iFit  = 0,
-                 Int_t iDrop = -1) {
+void pl_pull_trk(Int_t NSt = 8, Int_t iVar = 0, Int_t iFit = 0, Int_t iDrop = -1)
+{
   //  TCanvas *can = new TCanvas("can22","can22");
   //  can->Divide(2,2);
   TCanvas* can = new TCanvas("can", "can", 50, 0, 800, 800);
@@ -10,12 +8,14 @@ void pl_pull_trk(Int_t NSt   = 8,
     case 6:
     case 5:
     case 4: can->Divide(3, 3); break;
+    case 9: can->Divide(3, 4); break;
     case 18: can->Divide(4, 5); break;
     default: can->Divide(4, 4); ;
   }
   gPad->SetFillColor(0);
   gStyle->SetPalette(1);
   gStyle->SetOptStat(kTRUE);
+  gStyle->SetOptFit(kTRUE);
 
   gROOT->cd();
   gROOT->SetDirLevel(1);
@@ -57,12 +57,12 @@ void pl_pull_trk(Int_t NSt   = 8,
       gPad->SetLogy();
       gPad->SetGridx();
       if (iFit > 0) {
-        Double_t dFMean   = h1->GetMean();
-        Double_t dFLim    = 3.0 * h1->GetRMS();
-        Double_t dBinSize = h1->GetBinWidth(1);
-        dFLim             = TMath::Max(dFLim, 5. * dBinSize);
-        TFitResultPtr fRes =
-          h1->Fit("gaus", "S", "", dFMean - dFLim, dFMean + dFLim);
+        //Double_t dFMean   = h1->GetMean();
+        Double_t dFMean    = h1->GetBinCenter(h1->GetMaximumBin());
+        Double_t dFLim     = 2.0 * h1->GetRMS();
+        Double_t dBinSize  = h1->GetBinWidth(1);
+        dFLim              = TMath::Max(dFLim, 5. * dBinSize);
+        TFitResultPtr fRes = h1->Fit("gaus", "SQM", "", dFMean - dFLim, dFMean + dFLim);
         //cout << " fRes = "<< fRes <<endl;
         if (-1 == fRes) continue;
         if (iDrop == iSt) {  // drop station from deconvolution
@@ -78,7 +78,8 @@ void pl_pull_trk(Int_t NSt   = 8,
         //vSig[iIndSt]=TMath::Max(20.,vSig[iSt]);
         iIndSt++;
       }
-    } else {
+    }
+    else {
       cout << hname << " not found" << endl;
     }
   }
@@ -124,9 +125,8 @@ void pl_pull_trk(Int_t NSt   = 8,
   TMatrixD a(iIndSt, iIndSt);
   for (Int_t i = 0; i < iIndSt; i++)
     for (Int_t j = 0; j < iIndSt; j++) {
-      if (i == j) {
-        a[i][j] = 1;
-      } else {
+      if (i == j) { a[i][j] = 1; }
+      else {
         a[i][j] = 1. / val;
       }
     }
@@ -169,14 +169,8 @@ void pl_pull_trk(Int_t NSt   = 8,
   grr->Draw("APLE");
 
   for (Int_t i = 0; i < iIndSt; i++)
-    cout << Form(
-      "GMean %6.3f +/- %6.5f, GSig: %6.3f +/- %6.5f => ResC %d: %6.3f ",
-      vMean[i],
-      vMeanErr[i],
-      vSig[i],
-      vSigErr[i],
-      i,
-      vRes[i])
+    cout << Form("GMean %6.3f +/- %6.5f, GSig: %6.3f +/- %6.5f => ResC %d: %6.3f ", vMean[i], vMeanErr[i], vSig[i],
+                 vSigErr[i], i, vRes[i])
          << endl;
 
   cout << "Res-summary " << iVar << ": Nall, sigs = " << Nall;
diff --git a/reco/detectors/tof/CbmTofCalibrator.cxx b/reco/detectors/tof/CbmTofCalibrator.cxx
index f99d8dd4f5..2fb4dda54d 100644
--- a/reco/detectors/tof/CbmTofCalibrator.cxx
+++ b/reco/detectors/tof/CbmTofCalibrator.cxx
@@ -139,8 +139,8 @@ Bool_t CbmTofCalibrator::CreateCalHist()
   LOG(info) << "Define Calibrator histos for " << iNbDet << " detectors ";
 
   fhCalR0  = new TH1D("hCalR0", "Tracklet distance to nominal vertex; R_0 [cm]", 100, 0., 0.5);
-  fhCalDX0 = new TH1D("hCalDX0", "Tracklet distance to nominal vertex; #DeltaX_0 [cm]", 100, -0.5, 0.5);
-  fhCalDY0 = new TH1D("hCalDY0", "Tracklet distance to nominal vertex; #DeltaY_0 [cm]", 100, -0.5, 0.5);
+  fhCalDX0 = new TH1D("hCalDX0", "Tracklet distance to nominal vertex; #DeltaX_0 [cm]", 100, -1.5, 1.5);
+  fhCalDY0 = new TH1D("hCalDY0", "Tracklet distance to nominal vertex; #DeltaY_0 [cm]", 100, -1.5, 1.5);
 
   fhCalPosition.resize(iNbDet);
   fhCalPos.resize(iNbDet);
diff --git a/reco/detectors/tof/CbmTofEventClusterizer.cxx b/reco/detectors/tof/CbmTofEventClusterizer.cxx
index 1d7d110fe0..65fb304247 100644
--- a/reco/detectors/tof/CbmTofEventClusterizer.cxx
+++ b/reco/detectors/tof/CbmTofEventClusterizer.cxx
@@ -76,8 +76,9 @@ static Double_t dTLEvt        = 0.;
 static Int_t iNSpill          = 0;
 static Int_t iNbTs            = 0;
 
-const Double_t fdSpillDuration = 2.;   // in seconds
-const Double_t fdSpillBreak    = 0.9;  // in seconds
+const Double_t fdSpillDuration = 2.;    // in seconds
+const Double_t fdSpillBreak    = 0.9;   // in seconds
+const Double_t dTimeRes        = 0.08;  // in ns
 
 static Bool_t bAddBeamCounterSideDigi = kTRUE;
 static TRandom3* fRndm                = new TRandom3();
@@ -4712,6 +4713,7 @@ Bool_t CbmTofEventClusterizer::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc,
                                   vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                   //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                   Int_t(dLastTotS * 10.));  //channel -> Tot
+  pHit->SetTimeError(dTimeRes);
   // output hit
   new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
   if (fdMemoryTime > 0.) {  // memorize hit
@@ -5226,6 +5228,8 @@ Bool_t CbmTofEventClusterizer::BuildHits()
                                       vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                       //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                       Int_t(dWeightsSum * 10.));  //channel -> Tot
+                      pHit->SetTimeError(dTimeRes);
+
                       //0) ; //channel
                       // output hit
                       new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
@@ -5453,6 +5457,7 @@ Bool_t CbmTofEventClusterizer::BuildHits()
                                               vDigiIndRef.size(),  // number of linked digis =  2*CluSize
                                               //vPtsRef.size(), // flag  = number of TofPoints generating the cluster
                                               Int_t(dWeightsSum * 10.));  //channel -> Tot
+              pHit->SetTimeError(dTimeRes);
               //                0) ; //channel
               //                vDigiIndRef);
               // output hit
diff --git a/reco/detectors/tof/CbmTofFindTracks.cxx b/reco/detectors/tof/CbmTofFindTracks.cxx
index de064e29e8..2e2cc4197b 100644
--- a/reco/detectors/tof/CbmTofFindTracks.cxx
+++ b/reco/detectors/tof/CbmTofFindTracks.cxx
@@ -431,7 +431,7 @@ Bool_t CbmTofFindTracks::LoadCalParameter()
       new TH1F(Form("hPullT_Smt_Off"), Form("Tracklet PullT vs RpcInd ; RpcInd ; #DeltaT (ns)"), nSmt, 0, nSmt);
 
     // Initialize Parameter
-    if (fiCorMode == 3)  // hidden option, FIXME
+    if (fiCorMode <= 3)  // hidden option, FIXME
       for (Int_t iDet = 0; iDet < nSmt; iDet++) {
         std::map<Int_t, Int_t>::iterator it;
         //it = fMapRpcIdParInd.find(iDet);
@@ -650,14 +650,14 @@ Bool_t CbmTofFindTracks::WriteHistos()
             Double_t dLim      = 1.5 * hTOff1DY->GetRMS();
             TFitResultPtr fRes = hTOff1DY->Fit("gaus", "SQM", "", dMean - dLim, dMean + dLim);
             Int_t iFitStatus   = fRes;
-            if (iFitStatus == 0) {  // check validity of fit
-              dFMean = fRes->Parameter(1);
-              dVal += hTOff1D->GetBinContent(ix + 1);  //revert default correction
-              dVal -= dFMean;
-            }
-            LOG(info) << "Update hPullT_Smt_Off Ind " << ix << ": Old " << fhPullT_Smt_Off->GetBinContent(ix + 1)
-                      << ", Pull " << htmp1D->GetBinContent(ix + 1) << ", Dev@Peak " << hTOff1D->GetBinContent(ix + 1)
-                      << ", FitMean " << dFMean << " -> " << dVal;
+            //if (iFitStatus == 0) {  // check validity of fit
+            dFMean = fRes->Parameter(1);
+            dVal += hTOff1D->GetBinContent(ix + 1);  //revert default correction
+            dVal -= dFMean;
+            //}
+            LOG(info) << "Update hPullT_Smt_Off Ind " << ix << ", stat " << iFitStatus << ": Old "
+                      << fhPullT_Smt_Off->GetBinContent(ix + 1) << ", Pull " << htmp1D->GetBinContent(ix + 1)
+                      << ", Dev@Peak " << hTOff1D->GetBinContent(ix + 1) << ", FitMean " << dFMean << " -> " << dVal;
           }
           else {
             LOG(debug1) << "Update hPullT_Smt_Off " << ix << ": insufficient counts: " << hTOff1DY->GetEntries();
@@ -1253,8 +1253,8 @@ void CbmTofFindTracks::CreateHistograms()
   fhPullZ_Smt = new TH2F(Form("hPullZ_Smt"), Form("Tracklet ResiZ vs RpcInd ; RpcInd ; #DeltaZ (cm)"), nSmt, 0, nSmt,
                          100, -DZ0MAX, DZ0MAX);
 
-  fhTOff_Smt =
-    new TH2F(Form("hTOff_Smt"), Form("Tracklet TOff; RpcInd ; TOff (ns)"), nSmt, 0, nSmt, 501, -fT0MAX, fT0MAX);
+  fhTOff_Smt   = new TH2F(Form("hTOff_Smt"), Form("Tracklet TOff; RpcInd ; #DeltaTOff (ns)"), nSmt, 0, nSmt, 501,
+                        -fT0MAX * 20, fT0MAX * 20);
   fhTOff_HMul2 = new TH2F(Form("hTOff_HMul2"), Form("Tracklet TOff(HMul2); RpcInd ; TOff (ns)"), nSmt, 0, nSmt, 500,
                           -fT0MAX, fT0MAX);
 
-- 
GitLab