Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • le.koch/cbmroot
  • patrick.pfistner_AT_kit.edu/cbmroot
  • lena.rossel_AT_stud.uni-frankfurt.de/cbmroot
  • i.deppner/cbmroot
  • fweig/cbmroot
  • karpushkin_AT_inr.ru/cbmroot
  • v.akishina/cbmroot
  • rishat.sultanov_AT_cern.ch/cbmroot
  • l_fabe01_AT_uni-muenster.de/cbmroot
  • pwg-c2f/cbmroot
  • j.decuveland/cbmroot
  • a.toia/cbmroot
  • i.vassiliev/cbmroot
  • n.herrmann/cbmroot
  • o.lubynets/cbmroot
  • se.gorbunov/cbmroot
  • cornelius.riesen_AT_physik.uni-giessen.de/cbmroot
  • zhangqn17_AT_mails.tsinghua.edu.cn/cbmroot
  • bartosz.sobol/cbmroot
  • ajit.kumar/cbmroot
  • computing/cbmroot
  • a.agarwal_AT_vecc.gov.in/cbmroot
  • osingh/cbmroot
  • wielanek_AT_if.pw.edu.pl/cbmroot
  • malgorzata.karabowicz.stud_AT_pw.edu.pl/cbmroot
  • m.shiroya/cbmroot
  • s.roy/cbmroot
  • p.-a.loizeau/cbmroot
  • a.weber/cbmroot
  • ma.beyer/cbmroot
  • d.klein/cbmroot
  • d.smith/cbmroot
  • mvdsoft/cbmroot
  • d.spicker/cbmroot
  • y.h.leung/cbmroot
  • m.deveaux/cbmroot
  • mkunold/cbmroot
  • h.darwish/cbmroot
  • f_fido01_AT_uni-muenster.de/cbmroot
  • g.kozlov/cbmroot
  • d.emschermann/cbmroot
  • evgeny.lavrik/cbmroot
  • v.friese/cbmroot
  • f.uhlig/cbmroot
  • ebechtel_AT_ikf.uni-frankfurt.de/cbmroot
  • a.senger/cbmroot
  • praisig/cbmroot
  • s.lebedev/cbmroot
  • redelbach_AT_compeng.uni-frankfurt.de/cbmroot
  • p.subramani/cbmroot
  • a_meye37_AT_uni-muenster.de/cbmroot
  • om/cbmroot
  • o.golosov/cbmroot
  • l.chlad/cbmroot
  • a.bercuci/cbmroot
  • d.ramirez/cbmroot
  • v.singhal/cbmroot
  • h.schiller/cbmroot
  • apuntke/cbmroot
  • f.zorn/cbmroot
  • rubio_AT_physi.uni-heidelberg.de/cbmroot
  • p.chudoba/cbmroot
  • apuntke/mcbmroot
  • r.karabowicz/cbmroot
64 results
Show changes
Commits on Source (6)
Hide whitespace changes
Inline Side-by-side
macro/mcbm/mcbm_reco_event_tb_nh.C
View file @ aec14c0a
......@@ -13,7 +13,7 @@
// Binned tracker for track reconstruction
//
// V. Friese 11.06.2018
//
// S. Roy 11.01.2022 - added the Real event building and modified STS parAsic parameter
// --------------------------------------------------------------------------
void mcbm_reco_event_tb_nh(Int_t nEvents = 10, TString RunId = "test", TString InDir = "./data/",
......@@ -72,6 +72,20 @@ void mcbm_reco_event_tb_nh(Int_t nEvents = 10, TString RunId = "test", TString I
// setup->RemoveModule(ECbmModuleId::kSts);
// ------------------------------------------------------------------------
TString sEvBuildRaw = "Real";
// ----- Some global switches -----------------------------------------
// Bool_t eventBased = !sEvBuildRaw.IsNull();
Bool_t useMvd = setup->IsActive(ECbmModuleId::kMvd);
Bool_t useSts = setup->IsActive(ECbmModuleId::kSts);
Bool_t useRich = setup->IsActive(ECbmModuleId::kRich);
Bool_t useMuch = setup->IsActive(ECbmModuleId::kMuch);
Bool_t useTrd = setup->IsActive(ECbmModuleId::kTrd);
Bool_t useTof = setup->IsActive(ECbmModuleId::kTof);
Bool_t usePsd = setup->IsActive(ECbmModuleId::kPsd);
// ------------------------------------------------------------------------
// ----- Parameter files as input to the runtime database -------------
std::cout << std::endl;
......@@ -136,7 +150,7 @@ void mcbm_reco_event_tb_nh(Int_t nEvents = 10, TString RunId = "test", TString I
//mcManager->AddFile(rawFile);
//run->AddTask(mcManager);
// ------------------------------------------------------------------------
/*
CbmMcbm2018EventBuilder* eventBuilder = new CbmMcbm2018EventBuilder();
// eventBuilder->SetEventBuilderAlgo(EventBuilderAlgo::MaximumTimeGap);
// eventBuilder->SetMaximumTimeGap(100.);
......@@ -147,8 +161,60 @@ void mcbm_reco_event_tb_nh(Int_t nEvents = 10, TString RunId = "test", TString I
eventBuilder->SetTriggerMinNumberMuch(0);
eventBuilder->SetTriggerMinNumberTof(1);
eventBuilder->SetTriggerMinNumberRich(0);
eventBuilder->SetFillHistos(kTRUE);
if (timebased) run->AddTask(eventBuilder);
eventBuilder->SetFillHistos(kTRUE);*/
// ----- Raw event building from digis (the "Real" event builder) --------------------------------
if (sEvBuildRaw.EqualTo("Real", TString::ECaseCompare::kIgnoreCase)) {
CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
//Choose between NoOverlap, MergeOverlap, AllowOverlap
evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap);
// Remove detectors where digis not found
if (!useMvd) evBuildRaw->RemoveDetector(kRawEventBuilderDetMvd);
if (!useRich) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
if (!useMuch) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
if (!useTrd) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
if (!usePsd) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
if (!useTof) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
if (!useSts) {
std::cerr << "-E- " << myName << ": Sts must be present for raw event "
<< "building using ``Real2019'' option. Terminating macro." << std::endl;
return;
} // Set STS or Tof as reference detector
if (!useTof) evBuildRaw->SetReferenceDetector(kRawEventBuilderDetSts);
else
evBuildRaw->SetReferenceDetector(kRawEventBuilderDetTof);
// Use sliding window seed builder with STS
// evBuildRaw->SetReferenceDetector(kRawEventBuilderDetUndef);
// evBuildRaw->AddSeedTimeFillerToList(kRawEventBuilderDetSts);
// evBuildRaw->SetSlidingWindowSeedFinder(10, 40, 100);
// evBuildRaw->SetSeedFinderQa(true); // optional QA information for seed finder
evBuildRaw->SetTsParameters(0.0, 1.e7, 0.0); // Use CbmMuchDigi instead of CbmMuchBeamtimeDigi
evBuildRaw->ChangeMuchBeamtimeDigiFlag(kFALSE);
//Set event building parameters
if (!useTof) {
evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1);
evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1);
evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -10, 40);
}
else {
evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kTof, 1);
evBuildRaw->SetTriggerWindow(ECbmModuleId::kTof, -10, 40);
} //evBuildRaw->SetWriteHistosToFairSink(kFALSE);
//evBuildRaw->SetOutFilename("HistosEvtAllowOverlap_simulated.root");
run->AddTask(evBuildRaw);
std::cout << "-I- " << myName << ": Added task " << evBuildRaw->GetName() << std::endl;
} //? Real raw event building
// if (timebased) run->AddTask(eventBuilder);
// ----- Local reconstruction in MVD ----------------------------------
if (setup->IsActive(ECbmModuleId::kMvd)) {
......@@ -178,7 +244,7 @@ void mcbm_reco_event_tb_nh(Int_t nEvents = 10, TString RunId = "test", TString I
if (kFALSE && timebased) {
// ASIC params: #ADC channels, dyn. range, threshold, time resol., dead time,
// noise RMS, zero-threshold crossing rate
auto parAsic = new CbmStsParAsic(32, 75000., 3000., 5., 800., 1000., 3.9789e-3);
auto parAsic = new CbmStsParAsic(128, 32, 75000., 3000., 5., 800., 1000., 3.9789e-3);
// Module params: number of channels, number of channels per ASIC
auto parMod = new CbmStsParModule(2048, 128);
......
reco/KF/CbmKFParticleFinder.cxx
View file @ aec14c0a
......@@ -158,9 +158,17 @@ void CbmKFParticleFinder::Exec(Option_t* /*opt*/)
//calculate number of d-He4 candidates
int nCandidatesDHe4 = 0;
if (fPID)
for (int iTr = 0; iTr < nTracksEvent; iTr++)
if (TMath::Abs(fPID->GetPID()[iTr]) == 1000010029) nCandidatesDHe4++;
if (fPID) {
if ((int) fPID->GetPID().size() == nTracksEvent) {
for (int iTr = 0; iTr < nTracksEvent; iTr++) {
if (TMath::Abs(fPID->GetPID()[iTr]) == 1000010029) nCandidatesDHe4++;
}
}
else {
Error("CbmKFParticleFinder::Event", "PID task has a wrong number of tracks: %l of %l", fPID->GetPID().size(),
nTracksEvent);
}
}
vector<CbmStsTrack> vRTracks(nTracksEvent + nCandidatesDHe4);
vector<int> pdg(nTracksEvent + nCandidatesDHe4, -1);
......
reco/L1/CbmL1Performance.cxx
View file @ aec14c0a
......@@ -1025,6 +1025,22 @@ void CbmL1::TrackFitPerformance()
static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
static TH2F* pion_res_pt_fstt;
static TH2F* kaon_res_pt_fstt;
static TH2F* pton_res_pt_fstt;
static TH2F* pion_res_pt_vtxt;
static TH2F* kaon_res_pt_vtxt;
static TH2F* pton_res_pt_vtxt;
static TH2F* pion_res_p_fstt;
static TH2F* kaon_res_p_fstt;
static TH2F* pton_res_p_fstt;
static TH2F* pion_res_p_vtxt;
static TH2F* kaon_res_p_vtxt;
static TH2F* pton_res_p_vtxt;
static bool first_call = 1;
L1Fit fit;
......@@ -1043,6 +1059,22 @@ void CbmL1::TrackFitPerformance()
PRes2DPrim = new TH2F("PRes2DPrim", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
PRes2DSec = new TH2F("PRes2DSec", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
pion_res_pt_fstt = new TH2F("pion_res_pt_fstt", "", 100, 0, 10, 1000, -500, 500);
kaon_res_pt_fstt = new TH2F("kaon_res_pt_fstt", "", 100, 0, 10, 1000, -500, 500);
pton_res_pt_fstt = new TH2F("pton_res_pt_fstt", "", 100, 0, 10, 1000, -500, 500);
pion_res_pt_vtxt = new TH2F("pion_res_pt_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
kaon_res_pt_vtxt = new TH2F("kaon_res_pt_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
pton_res_pt_vtxt = new TH2F("pton_res_pt_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
pion_res_p_fstt = new TH2F("pion_res_p_fstt", "", 100, 0, 10, 1000, -500, 500);
kaon_res_p_fstt = new TH2F("kaon_res_p_fstt", "", 100, 0, 10, 1000, -500, 500);
pton_res_p_fstt = new TH2F("pton_res_p_fstt", "", 100, 0, 10, 1000, -500, 500);
pion_res_p_vtxt = new TH2F("pion_res_p_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
kaon_res_p_vtxt = new TH2F("kaon_res_p_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
pton_res_p_vtxt = new TH2F("pton_res_p_vtxt", "", 100, 0, 10, 1000, -5000, 5000);
struct {
const char* name;
const char* title;
......@@ -1151,8 +1183,14 @@ void CbmL1::TrackFitPerformance()
fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
L1Extrapolate(trPar, mcP.zIn, trPar.qp, fld);
h_fit[0]->Fill((trPar.x[0] - mcP.xIn) * 1.e4);
h_fit[1]->Fill((trPar.y[0] - mcP.yIn) * 1.e4);
double dx = trPar.x[0] - mcP.xIn;
double dy = trPar.y[0] - mcP.yIn;
double dt = sqrt(dx * dx + dy * dy);
// make dt distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
if (mc.ID % 2) dt = -dt;
double pt = sqrt(mcP.px * mcP.px + mcP.py * mcP.py);
h_fit[0]->Fill(dx * 1.e4);
h_fit[1]->Fill(dy * 1.e4);
h_fit[2]->Fill((trPar.tx[0] - mcP.pxIn / mcP.pzIn) * 1.e3);
h_fit[3]->Fill((trPar.ty[0] - mcP.pyIn / mcP.pzIn) * 1.e3);
h_fit[4]->Fill(fabs(1. / trPar.qp[0]) / mcP.p - 1);
......@@ -1160,9 +1198,26 @@ void CbmL1::TrackFitPerformance()
PRes2D->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
CbmL1MCTrack mcTrack = *(it->GetMCTracks()[0]);
if (mcTrack.IsPrimary()) PRes2DPrim->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
else
if (mcTrack.IsPrimary()) {
PRes2DPrim->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
if (abs(mcTrack.pdg == 211)) {
pion_res_p_fstt->Fill(mcP.p, dt * 1.e4);
pion_res_pt_fstt->Fill(pt, dt * 1.e4);
}
if (abs(mcTrack.pdg == 321)) {
kaon_res_p_fstt->Fill(mcP.p, dt * 1.e4);
kaon_res_pt_fstt->Fill(pt, dt * 1.e4);
}
if (abs(mcTrack.pdg == 2212)) {
pton_res_p_fstt->Fill(mcP.p, dt * 1.e4);
pton_res_pt_fstt->Fill(pt, dt * 1.e4);
}
}
else {
PRes2DSec->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
}
if (finite(trPar.C00[0]) && trPar.C00[0] > 0) h_fit[5]->Fill((trPar.x[0] - mcP.xIn) / sqrt(trPar.C00[0]));
if (finite(trPar.C11[0]) && trPar.C11[0] > 0) h_fit[6]->Fill((trPar.y[0] - mcP.yIn) / sqrt(trPar.C11[0]));
......@@ -1286,7 +1341,8 @@ void CbmL1::TrackFitPerformance()
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
L1TrackPar trPar(it->T, it->C);
if (mc.mother_ID != -1) { // secondary
// if (mc.mother_ID != -1) { // secondary
if (!mc.IsPrimary()) { // secondary
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
......@@ -1402,6 +1458,26 @@ void CbmL1::TrackFitPerformance()
}
if (mc.z != trPar.z[0]) continue;
double dx = trPar.x[0] - mc.x;
double dy = trPar.y[0] - mc.y;
double dt = sqrt(dx * dx + dy * dy);
// make dt distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
if (mc.ID % 2) dt = -dt;
double pt = sqrt(mc.px * mc.px + mc.py * mc.py);
if (abs(mc.pdg == 211)) {
pion_res_p_vtxt->Fill(mc.p, dt * 1.e4);
pion_res_pt_vtxt->Fill(pt, dt * 1.e4);
}
if (abs(mc.pdg == 321)) {
kaon_res_p_vtxt->Fill(mc.p, dt * 1.e4);
kaon_res_pt_vtxt->Fill(pt, dt * 1.e4);
}
if (abs(mc.pdg == 2212)) {
pton_res_p_vtxt->Fill(mc.p, dt * 1.e4);
pton_res_pt_vtxt->Fill(pt, dt * 1.e4);
}
// calculate pulls
h_fitPV[0]->Fill((mc.x - trPar.x[0]));
......@@ -1668,7 +1744,6 @@ void CbmL1::InputPerformance()
static TH1F *pullXmuch, *pullYmuch, *pullTmuch, *resXmuch, *resYmuch, *resTmuch /*, *pullZ*/;
static TH1F *pullXtrd, *pullYtrd, *pullTtrd, *resXtrd, *resYtrd, *resTtrd /*, *pullZ*/;
static TH1F *pullXtof, *pullYtof, *pullTtof, *resXtof, *resYtof, *resTtof /*, *pullZ*/;
static TH1F* trdMCPointsZ = nullptr;
static bool first_call = 1;
......@@ -1823,7 +1898,6 @@ void CbmL1::InputPerformance()
if (h.Det != 1) continue; // mvd hit
const CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(h.extIndex));
// int iMCPoint = -1;
CbmLink link;
CbmStsPoint* pt = 0;
......@@ -1847,19 +1921,15 @@ void CbmL1::InputPerformance()
}
}
if (stsHitMatch.GetNofLinks() > 0) {
Float_t bestWeight = 0.f;
for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
bestWeight = stsHitMatch.GetLink(iLink).GetWeight();
Int_t iFile = stsHitMatch.GetLink(iLink).GetFile();
Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry();
link = stsHitMatch.GetLink(iLink);
pt = (CbmStsPoint*) fStsPoints->Get(iFile, iEvent, stsHitMatch.GetLink(iLink).GetIndex());
link = stsHitMatch.GetLink(iLink);
bestWeight = link.GetWeight();
Int_t iFile = link.GetFile();
Int_t iEvent = link.GetEntry();
pt = (CbmStsPoint*) fStsPoints->Get(iFile, iEvent, link.GetIndex());
}
}
}
......@@ -1876,25 +1946,24 @@ void CbmL1::InputPerformance()
sh->Position(hitPos);
sh->PositionError(hitErr);
// pt->Position(mcPos); // this is wrong!
mcPos.SetX(pt->GetX(hitPos.Z()));
mcPos.SetY(pt->GetY(hitPos.Z()));
double t = (hitPos.Z() - pt->GetZIn()) / (pt->GetZOut() - pt->GetZIn());
mcPos.SetX(pt->GetXIn() + t * (pt->GetXOut() - pt->GetXIn()));
mcPos.SetY(pt->GetYIn() + t * (pt->GetYOut() - pt->GetYIn()));
mcPos.SetZ(hitPos.Z());
#if 0 // standard errors
if (hitErr.X() != 0) pullXsts->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
if (hitErr.Y() != 0) pullYsts->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
#elif 1 // qa errors
if (hitErr.X() != 0) pullXsts->Fill((hitPos.X() - mcPos.X()) / sh->GetDx()); // qa errors
if (hitErr.Y() != 0) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
#else // errors used in TF
if (hitErr.X() != 0)
if (0) { // standard errors
if (hitErr.X() > 1.e-8) pullXsts->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
if (hitErr.Y() > 1.e-8) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
}
else if (1) { // qa errors
if (sh->GetDx() > 1.e-8) pullXsts->Fill((hitPos.X() - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
}
else { // errors used in TF
pullXsts->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
if (hitErr.Y() != 0)
pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
#endif
}
resXsts->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
resYsts->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
......@@ -1995,20 +2064,19 @@ void CbmL1::InputPerformance()
mcPos.SetY(0.5 * (pt->GetYIn() + pt->GetYOut()));
mcPos.SetZ(hitPos.Z());
#if 0 // standard errors
if (hitErr.X() != 0) pullXmuch->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
if (hitErr.Y() != 0) pullYmuch->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
#elif 1 // qa errors
if (hitErr.X() != 0) pullXmuch->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
if (hitErr.Y() != 0) pullYmuch->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTmuch->Fill((h.t - mcTime) / sh->GetTimeError());
#else // errors used in TF
if (hitErr.X() != 0)
if (0) { // standard errors
if (hitErr.X() > 1.e-8) pullXmuch->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
if (hitErr.Y() > 1.e-8) pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
}
else if (1) { // qa errors
if (sh->GetDx() > 1.e-8) pullXmuch->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYmuch->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTmuch->Fill((h.t - mcTime) / sh->GetTimeError());
}
else { // errors used in TF
pullXmuch->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
if (hitErr.Y() != 0)
pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
#endif
}
resXmuch->Fill((h.x - mcPos.X()) * 10 * 1000);
resYmuch->Fill((h.y - mcPos.Y()) * 10 * 1000);
......@@ -2027,12 +2095,13 @@ void CbmL1::InputPerformance()
if (hm->GetNofLinks() == 0) continue;
if (hm->GetNofLinks() != 1) continue; //SG!!
if (hm->GetNofLinks() != 1) continue; // only check single-linked hits
Float_t bestWeight = 0.f;
Float_t totalWeight = 0.f;
int iMCPoint = -1;
CbmLink link;
for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
totalWeight += hm->GetLink(iLink).GetWeight();
if (hm->GetLink(iLink).GetWeight() > bestWeight) {
......@@ -2062,18 +2131,19 @@ void CbmL1::InputPerformance()
mcPos.SetY((pt->GetYIn() + pt->GetYOut()) / 2.);
mcPos.SetZ(hitPos.Z());
#if 0 // standard errors
if (hitErr.X() != 0) pullXtrd->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
if (hitErr.Y() != 0) pullYtrd->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
#elif 1 // qa errors
if (hitErr.X() != 0) pullXtrd->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
if (hitErr.Y() != 0) pullYtrd->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtrd->Fill((h.t - mcTime) / sh->GetTimeError());
#else // errors used in TF
if (hitErr.X() != 0) pullXtrd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
if (hitErr.Y() != 0) pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
#endif
if (0) { // standard errors
if (hitErr.X() > 1.e-8) pullXtrd->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
if (hitErr.Y() > 1.e-8) pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
}
else if (1) { // qa errors
if (sh->GetDx() > 1.e-8) pullXtrd->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYtrd->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtrd->Fill((h.t - mcTime) / sh->GetTimeError());
}
else { // errors used in TF
pullXtrd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
}
resXtrd->Fill((h.x - mcPos.X()) * 10 * 1000);
resYtrd->Fill((h.y - mcPos.Y()) * 10 * 1000);
......@@ -2129,18 +2199,19 @@ void CbmL1::InputPerformance()
mcPos.SetY((pt->GetY()));
mcPos.SetZ(hitPos.Z());
#if 0 // standard errors
if (hitErr.X() != 0) pullXmuch->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
if (hitErr.Y() != 0) pullYmuch->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
#elif 1 // qa errors
if (hitErr.X() != 0) pullXtof->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
if (hitErr.Y() != 0) pullYtof->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
#else // errors used in TF
if (hitErr.X() != 0) pullXtof->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
if (hitErr.Y() != 0) pullYtof->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
#endif
if (0) { // standard errors
if (hitErr.X() > 1.e-8) pullXmuch->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
if (hitErr.Y() > 1.e-8) pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
}
else if (1) { // qa errors
if (sh->GetDx() > 1.e-8) pullXtof->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYtof->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
}
else { // errors used in TF
pullXtof->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
pullYtof->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
}
resXtof->Fill((h.x - mcPos.X()) * 10 * 1000);
resYtof->Fill((h.y - mcPos.Y()) * 10 * 1000);
......
reco/detectors/sts/CbmRecoStsPixel.cxx
View file @ aec14c0a
......@@ -72,7 +72,7 @@ InitStatus CbmRecoStsPixel::Init()
{
// --- Something for the screen
std::cout << std::endl;
LOG(info) << "==========================================================";
LOG(info) << GetName() << ": Initialising ";
......@@ -259,31 +259,28 @@ void CbmRecoStsPixel::ProcessData(CbmEvent* event)
break;
}
double dx = pt->GetXOut() - pt->GetXIn();
double dy = pt->GetYOut() - pt->GetYIn();
double dz = pt->GetZOut() - pt->GetZIn();
if (fabs(dz) > 1.e-2) {
hit->SetX(pt->GetXIn() + dx * (staZ - pt->GetZIn()) / dz);
hit->SetY(pt->GetYIn() + dy * (staZ - pt->GetZIn()) / dz);
hit->SetZ(staZ);
}
else {
hit->SetX(pt->GetXIn());
hit->SetY(pt->GetYIn());
hit->SetZ(pt->GetZIn());
}
hit->SetX(0.5 * (pt->GetXIn() + pt->GetXOut()));
hit->SetY(0.5 * (pt->GetYIn() + pt->GetYOut()));
hit->SetZ(0.5 * (pt->GetZIn() + pt->GetZOut()));
Double_t pitchX = station->GetSensorPitch(0);
Double_t pitchY = station->GetSensorPitch(1);
Double_t resX = station->GetSensorPitch(0) / sqrt(12.);
Double_t resY = station->GetSensorPitch(1) / sqrt(12.);
assert(pitchX > 1.e-5);
assert(pitchY > 1.e-5);
assert(resX > 1.e-5);
assert(resY > 1.e-5);
hit->SetX((floor(hit->GetX() / pitchX) + 0.5) * pitchX);
hit->SetY((floor(hit->GetY() / pitchY) + 0.5) * pitchY);
auto gaus = []() {
double x = 5;
while (fabs(x) > 3.5) {
x = gRandom->Gaus(0., 1.);
}
return x;
};
hit->SetDx(pitchX / sqrt(12.));
hit->SetDy(pitchY / sqrt(12.));
hit->SetX(hit->GetX() + resX * gaus());
hit->SetY(hit->GetY() + resY * gaus());
hit->SetDx(resX);
hit->SetDy(resY);
hit->SetDxy(0.);
hit->SetDu(hit->GetDx());
hit->SetDv(hit->GetDy());
......
sim/detectors/sts/CbmStsDigitizePixel.cxx
View file @ aec14c0a
......@@ -74,11 +74,17 @@ using namespace CbmSts;
ClassImp(CbmStsDigitizePixel);
// ----- Standard constructor ------------------------------------------
CbmStsDigitizePixel::CbmStsDigitizePixel(Double_t pitchXcm, Double_t pitchYcm, Double_t resolutionTns)
CbmStsDigitizePixel::CbmStsDigitizePixel() : CbmDigitize<CbmStsDigi>("StsDigitizePixel") {}
// -------------------------------------------------------------------------
// ----- Standard constructor ------------------------------------------
CbmStsDigitizePixel::CbmStsDigitizePixel(Double_t resolutionXcm, Double_t resolutionYcm, Double_t resolutionTns,
int nPixelStations)
: CbmDigitize<CbmStsDigi>("StsDigitizePixel")
, fPitchXcm(pitchXcm)
, fPitchYcm(pitchYcm)
, fResolutionTns(resolutionTns)
, fPixelNstations(nPixelStations)
, fPixelResolutionXcm(resolutionXcm)
, fPixelResolutionYcm(resolutionYcm)
, fPixelResolutionTns(resolutionTns)
{
}
// -------------------------------------------------------------------------
......@@ -175,7 +181,12 @@ void CbmStsDigitizePixel::Exec(Option_t* /*opt*/)
UInt_t address = static_cast<UInt_t>(point->GetDetectorID());
UShort_t channel = 0;
double timef = fCurrentEventTime + point->GetTime();
timef += gRandom->Gaus(0, fResolutionTns);
if (fSetup->GetStationNumber(address) < fPixelNstations) { timef += gRandom->Gaus(0, fPixelResolutionTns); }
else {
timef += gRandom->Gaus(0, fStripResolutionTns);
}
Long64_t time = Long64_t(round(timef));
if (time < 0) { time = 0; }
UShort_t adc = 30;
......@@ -195,7 +206,7 @@ void CbmStsDigitizePixel::Exec(Option_t* /*opt*/)
}
timer.Stop();
if (fCreateMatches) { std::cout << "create matches" << std::endl; }
// --- Event log
LOG(info) << left << setw(15) << GetName() << "[" << fixed << setprecision(3) << timer.RealTime() << " s]"
<< " Points processed " << fPoints->GetEntriesFast();
......@@ -231,28 +242,40 @@ void CbmStsDigitizePixel::InitParams()
UInt_t nAsicChannels = 128; // Number of readout channels per ASIC
// --- ASIC parameters
UShort_t nAdc = 32; // Number of ADC channels (5 bit)
Double_t dynRange = 75000.; // Dynamic range [e]
Double_t threshold = 3000.; // Threshold [e]
Double_t timeResol = fResolutionTns; // Time resolution [ns]
Double_t deadTime = 800.; // Channel dead time [ns]
Double_t noiseRms = 1000.; // RMS of noise [e]
Double_t znr = 3.9789e-3; // Zero-crossing noise rate [1/ns]
CbmStsParAsic userParAsic(nAsicChannels, nAdc, dynRange, threshold, timeResol, deadTime, noiseRms, znr);
CbmStsParModule userParModule(nChannels, nAsicChannels);
userParModule.SetAllAsics(userParAsic);
UShort_t nAdc = 32; // Number of ADC channels (5 bit)
Double_t dynRange = 75000.; // Dynamic range [e]
Double_t threshold = 3000.; // Threshold [e]
Double_t deadTime = 800.; // Channel dead time [ns]
Double_t noiseRms = 1000.; // RMS of noise [e]
Double_t znr = 3.9789e-3; // Zero-crossing noise rate [1/ns]
CbmStsParAsic userParAsicStrip(nAsicChannels, nAdc, dynRange, threshold, fStripResolutionTns, deadTime, noiseRms,
znr);
CbmStsParModule userParModuleStrip(nChannels, nAsicChannels);
userParModuleStrip.SetAllAsics(userParAsicStrip);
CbmStsParAsic userParAsicPixel(nAsicChannels, nAdc, dynRange, threshold, fPixelResolutionTns, deadTime, noiseRms,
znr);
CbmStsParModule userParModulePixel(nChannels, nAsicChannels);
userParModulePixel.SetAllAsics(userParAsicPixel);
assert(fParSetModule);
fParSetModule->clear();
for (Int_t iModule = 0; iModule < fSetup->GetNofModules(); iModule++) {
UInt_t address = fSetup->GetModule(iModule)->GetAddress();
fParSetModule->SetParModule(address, userParModule);
if (fSetup->GetStationNumber(address) < fPixelNstations) {
fParSetModule->SetParModule(address, userParModulePixel);
}
else {
fParSetModule->SetParModule(address, userParModuleStrip);
}
}
fParSetModule->setChanged();
fParSetModule->setInputVersion(-2, 1);
LOG(info) << GetName() << "--- Using global ASIC parameters: \n " << userParAsic.ToString();
LOG(info) << GetName() << "--- Using global ASIC parameters for Pixels: \n " << userParAsicPixel.ToString();
LOG(info) << GetName() << "--- Using global ASIC parameters for Strips: \n " << userParAsicPixel.ToString();
LOG(info) << GetName() << "--- Module parameters: " << fParSetModule->ToString();
}
......@@ -287,10 +310,17 @@ void CbmStsDigitizePixel::InitParams()
// sensor parameters used by the tracker
par.SetPar(6, fPitchXcm); // Strip pitch front side
par.SetPar(7, fPitchYcm); // Strip pitch back side
par.SetPar(8, 0.); // Stereo angle front side [deg]
par.SetPar(9, 90); // Stereo angle back side [deg]
if (fSetup->GetStationNumber(address) < fPixelNstations) {
par.SetPar(6, fPixelResolutionXcm * sqrt(12.)); // Strip pitch front side
par.SetPar(7, fPixelResolutionYcm * sqrt(12.)); // Strip pitch back side
}
else {
par.SetPar(6, fStripResolutionXcm * sqrt(12.)); // Strip pitch front side
par.SetPar(7, fStripResolutionYcm * sqrt(12.)); // Strip pitch back side
}
par.SetPar(8, 0.); // Stereo angle front side [deg]
par.SetPar(9, 90); // Stereo angle back side [deg]
fParSetSensor->SetParSensor(address, par);
}
......
sim/detectors/sts/CbmStsDigitizePixel.h
View file @ aec14c0a
......@@ -54,8 +54,10 @@ class CbmStsDigitizePixel : public CbmDigitize<CbmStsDigi> {
public:
/** Constructor **/
CbmStsDigitizePixel(Double_t pitchXcm = 0.0010, Double_t pitchYcm = 0.0010, Double_t resolutionTns = 5.);
CbmStsDigitizePixel();
/** Constructor **/
CbmStsDigitizePixel(Double_t resolutionXcm, Double_t resolutionYcm, Double_t resolutionTns, int nPixelStations);
/** Destructor **/
virtual ~CbmStsDigitizePixel();
......@@ -86,14 +88,28 @@ public:
/** Initialisation **/
virtual InitStatus Init();
/// set pitch X [cm]
void SetPitchX(double pitchXcm) { fPitchXcm = pitchXcm; }
/// set number of pixel stations
void SetPixelNstations(int n) { fPixelNstations = n; }
/// set pixel resolution X [cm]
void SetPixelResolutionX(double resXcm) { fPixelResolutionXcm = resXcm; }
/// set pixel resolution Y [cm]
void SetPixelResolutionY(double resYcm) { fPixelResolutionYcm = resYcm; }
/// set pixel resolution Time [ns]
void SetPixelResolutionTime(double resolutionTns) { fPixelResolutionTns = resolutionTns; }
/// set pitch Y [cm]
void SetPitchY(double pitchYcm) { fPitchYcm = pitchYcm; }
/// set strip resolution X [cm]
void SetStripResolutionX(double resXcm) { fStripResolutionXcm = resXcm; }
/// set strip resolution Y [cm]
void SetStripResolutionY(double resYcm) { fStripResolutionYcm = resYcm; }
/// set strip resolution Time [ns]
void SetStripResolutionTime(double resolutionTns) { fStripResolutionTns = resolutionTns; }
/// set resolution Time [ns]
void SetResolutionTime(double resolutionTns) { fResolutionTns = resolutionTns; }
private:
Bool_t fIsInitialised; ///< kTRUE if Init() was called
......@@ -110,9 +126,17 @@ private:
// data members
Double_t fPitchXcm {0.0058}; // pitch in X [cm]
Double_t fPitchYcm {0.0058}; // pitch in Y [cm]
Double_t fResolutionTns {5.}; // resolution in time [ns]
Int_t fPixelNstations {3}; // number of pixel stations.
// First fPixelNStations stations of STS will be replaced by pixels.
Double_t fPixelResolutionXcm {0.0010 / sqrt(12.)}; // pixel resolution in X [cm]
Double_t fPixelResolutionYcm {0.0010 / sqrt(12.)}; // pixel resolution in Y [cm]
Double_t fPixelResolutionTns {2.6}; // pixel resolution in time [ns]
Double_t fStripResolutionXcm {0.0010}; // strip resolution in X [cm]
Double_t fStripResolutionYcm {0.0100}; // strip resolution in Y [cm]
Double_t fStripResolutionTns {2.6}; // strip resolution in time [ns]
/** @brief Initialise the parameters **/
void InitParams();
......