From e1719e4c6c191584771b69f8bd7e972e0cecd229 Mon Sep 17 00:00:00 2001
From: Valentina <v.akishina@gsi.de>
Date: Wed, 23 Feb 2022 12:55:15 +0100
Subject: [PATCH] L1: add option for mc-point-based reconstruction
---
reco/L1/CbmL1.cxx | 32 ++-
reco/L1/CbmL1.h | 14 ++
reco/L1/CbmL1MCPoint.h | 1 +
reco/L1/CbmL1ReadEvent.cxx | 387 +++++++++++++++++++++----------------
reco/L1/L1Algo/L1Algo.cxx | 2 +
reco/L1/L1Algo/L1Station.h | 2 +-
6 files changed, 261 insertions(+), 177 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 96123a2d7c..ec688711b4 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -195,7 +195,7 @@ InitStatus CbmL1::Init()
fUseTOF = 0;
if (fTrackingMode == L1Algo::TrackingMode::kMcbm) {
- fUseMUCH = 0;
+ fUseMUCH = 1;
fUseTRD = 1;
fUseTOF = 1;
}
@@ -433,6 +433,8 @@ InitStatus CbmL1::Init()
if (NULL == fChannelInfo) break;
float z = fChannelInfo->GetZ();
float x = fChannelInfo->GetX();
+ float y = fChannelInfo->GetY();
+
if (fMissingHits) {
if ((x > 20) && (z > 270) && (station == 1)) station = 2;
if (z > 400) continue;
@@ -534,11 +536,12 @@ InitStatus CbmL1::Init()
geo.push_back(t.RadLength);
fscal f_phi = 0, f_sigma = mvdStationPar->GetXRes(ist) / 10000, b_phi = 3.14159265358 / 2.,
- b_sigma = mvdStationPar->GetYRes(ist) / 10000;
+ b_sigma = mvdStationPar->GetYRes(ist) / 10000, dt = 1000;
geo.push_back(f_phi);
geo.push_back(f_sigma);
geo.push_back(b_phi);
geo.push_back(b_sigma);
+ geo.push_back(dt);
z = t.z;
Xmax = Ymax = t.R;
@@ -564,6 +567,7 @@ InitStatus CbmL1::Init()
b_phi += station->GetSensorStereoAngle(1) * Pi / 180.;
f_sigma = station->GetSensorPitch(0) / TMath::Sqrt(12);
b_sigma = f_sigma;
+ fscal dt = 5;
//f_sigma *= cos(f_phi); // TODO: think about this
//b_sigma *= cos(b_phi);
@@ -571,6 +575,7 @@ InitStatus CbmL1::Init()
geo.push_back(f_sigma);
geo.push_back(b_phi);
geo.push_back(b_sigma);
+ geo.push_back(dt);
z = station->GetZ();
Xmax = station->GetXmax();
@@ -597,15 +602,17 @@ InitStatus CbmL1::Init()
geo.push_back(2);
geo.push_back(z);
geo.push_back(layer->GetDz());
- geo.push_back(0);
+ geo.push_back(10);
geo.push_back(100); //station->GetRmax()
geo.push_back(0);
- fscal f_phi = 0, f_sigma = 0.1, b_phi = 3.14159265358 / 2., b_sigma = 0.1;
+ // fscal f_phi = 0, f_sigma = 0.1, b_phi = 3.14159265358 / 2., b_sigma = 0.1;
+ fscal f_phi = 0, f_sigma = 0.35, b_phi = 3.14159265358 / 2., b_sigma = 0.35, dt = 3.9;
geo.push_back(f_phi);
geo.push_back(f_sigma);
geo.push_back(b_phi);
geo.push_back(b_sigma);
+ geo.push_back(dt);
Xmax = 100; //station->GetRmax();
@@ -644,15 +651,17 @@ InitStatus CbmL1::Init()
geo.push_back(stationZ);
geo.push_back(2 * module->GetSizeZ());
- geo.push_back(0);
+ geo.push_back(10);
geo.push_back(2 * module->GetSizeX());
geo.push_back(10);
- fscal f_phi = 0, f_sigma = 1., b_phi = 3.14159265358 / 2., b_sigma = 1.;
+ // fscal f_phi = 0, f_sigma = 1., b_phi = 3.14159265358 / 2., b_sigma = 1.;
+ fscal f_phi = 0, f_sigma = 0.15, b_phi = 3.14159265358 / 2., b_sigma = 0.15, dt = 10;
geo.push_back(f_phi);
geo.push_back(f_sigma);
geo.push_back(b_phi);
geo.push_back(b_sigma);
+ geo.push_back(dt);
Xmax = module->GetSizeX();
Ymax = module->GetSizeY();
LOG(info) << "L1: Trd station " << num << " at z " << stationZ << endl;
@@ -672,11 +681,13 @@ InitStatus CbmL1::Init()
geo.push_back(150); /// TODO: add Tof max radius
geo.push_back(10);
- fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., b_sigma = 1 / 10000;
+ // fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., b_sigma = 1 / 10000;
+ fscal f_phi = 0, f_sigma = 0.42, b_phi = 3.14159265358 / 2., b_sigma = 0.23, dt = 0.075;
geo.push_back(f_phi);
geo.push_back(f_sigma);
geo.push_back(b_phi);
geo.push_back(b_sigma);
+ geo.push_back(dt);
Xmax = Ymax = 20;
LOG(info) << "L1: Tof station " << num << " at z " << z << endl;
}
@@ -875,7 +886,7 @@ InitStatus CbmL1::Init()
// MVD // TODO: to be exchanged with specific flags (timeInfo, fieldInfo etc.) (S.Zh.)
stationInfo.SetTimeInfo(1);
stationInfo.SetZ(layer->GetZ());
- stationInfo.SetMaterial(layer->GetDz(), 0); // TODO: Why rad len is 0????? (S.Zh.)
+ stationInfo.SetMaterial(layer->GetDz(), 10); // TODO: Why rad len is 0????? (S.Zh.)
stationInfo.SetXmax(100.);
stationInfo.SetYmax(100.);
stationInfo.SetRmin(0.);
@@ -1073,10 +1084,11 @@ InitStatus CbmL1::Init()
const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
algo->fRadThick[iSta].SetBins(NBins, RMax);
algo->fRadThick[iSta].table.resize(NBins);
+
+ float hole = 0.15;
for (int iB = 0; iB < NBins; iB++) {
algo->fRadThick[iSta].table[iB].resize(NBins);
- float hole = 0.15;
for (int iB2 = 0; iB2 < NBins; iB2++) {
algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
// Correction for holes in material map
@@ -1137,10 +1149,10 @@ InitStatus CbmL1::Init()
const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
algo->fRadThick[iSta].SetBins(NBins, RMax);
algo->fRadThick[iSta].table.resize(NBins);
+ float hole = 0.15;
for (int iB = 0; iB < NBins; iB++) {
algo->fRadThick[iSta].table[iB].resize(NBins);
- float hole = 0.15;
for (int iB2 = 0; iB2 < NBins; iB2++) {
algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
// Correction for holes in material map
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 8be4660f0e..75bb1aac03 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -149,6 +149,14 @@ public:
/// Sets a vector of detectors used in tracking
void SetActiveTrackingDetectorIDs(const std::set<L1DetectorID>& init) { fActiveTrackingDetectorIDs = init; }
+ void SetUseMcHit(int StsUseMcHit = 0, int MuchUseMcHit = 0, int TrdUseMcHit = 0, int TofUseMcHit = 0)
+ {
+ fStsUseMcHit = StsUseMcHit;
+ fMuchUseMcHit = MuchUseMcHit;
+ fTrdUseMcHit = TrdUseMcHit;
+ fTofUseMcHit = TofUseMcHit;
+ }
+
void SetStsMaterialBudgetFileName(TString fileName) { fStsMatBudgetFileName = fileName; }
void SetMvdMaterialBudgetFileName(TString fileName) { fMvdMatBudgetFileName = fileName; }
@@ -278,6 +286,11 @@ private:
Double_t fMomentumCutOff {0.1}; // currently not used
Bool_t fGhostSuppression {true}; // currently not used
+ Int_t fStsUseMcHit {-1}; // if STS data should be processed
+ Int_t fMuchUseMcHit {-1}; // if Much data should be processed
+ Int_t fTrdUseMcHit {-1}; // if Trd data should be processed
+ Int_t fTofUseMcHit {-1}; // if Tof data should be processed
+
Bool_t fUseMVD {false}; // if Mvd data should be processed
Bool_t fUseMUCH {false}; // if Much data should be processed
Bool_t fUseTRD {false}; // if Trd data should be processed
@@ -326,6 +339,7 @@ private:
TClonesArray* fTofHits {nullptr}; // CbmMatches array
CbmTofDigiPar* fDigiPar {nullptr};
CbmTofDigiBdfPar* fTofDigiBdfPar {nullptr};
+ vector<vector<int>> TofPointToTrack;
TFile* fPerfFile {nullptr};
TDirectory* fHistoDir {nullptr};
diff --git a/reco/L1/CbmL1MCPoint.h b/reco/L1/CbmL1MCPoint.h
index bb2bf315b8..660d1635bb 100644
--- a/reco/L1/CbmL1MCPoint.h
+++ b/reco/L1/CbmL1MCPoint.h
@@ -66,6 +66,7 @@ struct CbmL1MCPoint {
int pointId = -1;
int file = -1;
int event = -1;
+ int trackId = -1;
L1Vector<int> hitIds {"CbmL1MCPoint::hitIds"}; // indices of CbmL1Hits in L1->vStsHits array
};
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 4c8395355a..a32954f561 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -74,17 +74,71 @@ struct TmpHit { // used for sort Hits before writing in the normal arrays
int ExtIndex; // index of hit in the TClonesArray array ( negative for MVD )
double u_front, u_back; // positions of strips
double x, y, z; // position of hit
+ double xmc, ymc, p; // mc-position of hit, momentum
+ double tx, ty; // mc-slope of mc point
double dx, dy, dxy;
double du, dv;
int iMC; // index of MCPoint in the vMCPoints array
double time = 0., dt = 1.e10;
int Det;
int id;
+ int track;
static bool Compare(const TmpHit& a, const TmpHit& b)
{
return (a.iStation < b.iStation) || ((a.iStation == b.iStation) && (a.y < b.y));
}
+
+ void CreateHitFromPoint(const CbmL1MCPoint& point, int det, int nTmpHits, int nStripF, int ip, int& NStrips,
+ L1Station& st)
+ {
+ ExtIndex = 0;
+ Det = det;
+ id = nTmpHits; //tmpHits.size();
+ iStation = point.iStation;
+
+ dt = st.dt[0];
+ time = point.time; // + gRandom->Gaus(0, dt);
+
+ iStripF = nStripF + ip; //firstDetStrip + ip;
+ iStripB = iStripF;
+ if (NStrips <= iStripF) { NStrips = iStripF + 1; }
+
+ float f_sigma = 0.1; // sqrt(st.frontInfo.sigma2[0]);
+ //1.;
+ float b_sigma = 0.1; // sqrt(st.backInfo.sigma2[0]);//1.;
+
+ dx = f_sigma;
+ dy = b_sigma;
+ dxy = 0;
+ du = f_sigma;
+ dv = b_sigma;
+
+ x = point.x + gRandom->Gaus(0, dx);
+ y = point.y + gRandom->Gaus(0, dy);
+ z = point.z;
+
+ xmc = point.x;
+ ymc = point.y;
+
+ track = point.ID;
+
+ p = point.p;
+
+ u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
+ u_back = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
+ iMC = ip;
+ }
+
+ void SetHitFromPoint(const CbmL1MCPoint& point, L1Station& st)
+ {
+ x = 0.5 * (point.xIn + point.xOut) + gRandom->Gaus(0, dx);
+ y = 0.5 * (point.yIn + point.yOut) + gRandom->Gaus(0, dy);
+ z = 0.5 * (point.zIn + point.zOut);
+ time = point.time + gRandom->Gaus(0, dt);
+ u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
+ u_back = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
+ }
};
/// Repack data from Clones Arrays to L1 arrays
@@ -153,9 +207,12 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
// get STS hits
int nStsHits = 0;
- L1Vector<CbmLink*> ToFPointsMatch("CbmL1ReadEvent::ToFPointsMatch");
+ int firstTrdPoint = 0;
+ int firstStsPoint = 0;
+ int firstMuchPoint = 0;
+ int firstTofPoint = 0;
- int firstTrdPoint = 0;
+ // L1Vector<CbmLink*> ToFPointsMatch("CbmL1ReadEvent::ToFPointsMatch");
if (fPerformance) {
Fill_vMCTracks();
@@ -202,10 +259,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
}
if (fVerbose > 2) { LOG(info) << "CbmL1ReadEvent: max deviation of Mvd points " << maxDeviation; }
// ensure that the nominal station positions are not far from the sensors
- assert(fabs(maxDeviation) < 1.);
+ // assert(fabs(maxDeviation) < 1.);
}
-
+ firstStsPoint = vMCPoints.size();
if (fStsPoints) {
Int_t nMC = fStsPoints->Size(iFile, iEvent);
double maxDeviation = 0;
@@ -239,9 +296,10 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
}
if (fVerbose > 2) { LOG(info) << "CbmL1ReadEvent: max deviation of Sts points " << maxDeviation; }
// ensure that the nominal station positions are not far from the sensors
- assert(fabs(maxDeviation) < 1.);
+ // assert(fabs(maxDeviation) < 1.);
}
+ firstMuchPoint = vMCPoints.size();
if (fMuchPoints) {
Int_t nMC = fMuchPoints->Size(iFile, iEvent);
for (Int_t iMC = 0; iMC < nMC; iMC++) {
@@ -292,67 +350,37 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
}
}
- ToFPointsMatch.clear();
+ //ToFPointsMatch.clear();
+ firstTofPoint = vMCPoints.size();
if (fTofPoints) {
- // TOF data arrays also contain fake beam-counter data. Select the real TOF points here.
-
- for (int j = 0; j < fTofHits->GetEntriesFast(); j++) {
-
- CbmLink* link = 0;
-
- CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j));
-
- if (0x00202806 == mh->GetAddress() || 0x00002806 == mh->GetAddress()) {
- ToFPointsMatch.push_back(link);
- continue;
- }
-
- CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(j));
-
- if (matchHitMatch->GetNofLinks() > 0) {
-
- link = (CbmLink*) &matchHitMatch->GetLink(0);
-
- CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(link->GetFile(), link->GetEntry(), link->GetIndex());
-
- for (int iLink = 1; iLink < matchHitMatch->GetNofLinks(); iLink++) {
-
- CbmLink* link1 = (CbmLink*) &(matchHitMatch->GetLink(iLink));
-
- CbmTofPoint* pt_cur =
- (CbmTofPoint*) fTofPoints->Get(link1->GetFile(), link1->GetEntry(), link1->GetIndex());
- TVector3 pos_cur, pos_old, pos_hit;
+ vector<float> TofPointToTrackdZ[NTOFStation];
+ TofPointToTrack.resize(NTOFStation);
- pt_cur->Position(pos_cur);
- pt->Position(pos_old);
- mh->Position(pos_hit);
+ for (Int_t i = 0; i < NTOFStation; i++) {
- if (fabs(pos_cur.Z() - pos_hit.Z()) < fabs(pos_old.Z() - pos_hit.Z())) {
- pt = pt_cur;
- link = link1;
- }
- }
-
- ToFPointsMatch.push_back(link);
- } // for j
- } // if listTrdHits
+ TofPointToTrack[i].resize(vMCTracks.size(), 0);
+ TofPointToTrackdZ[i].resize(vMCTracks.size());
+ }
+ for (int iSt = 0; iSt < NTOFStation; iSt++)
+ for (unsigned int i = 0; i < TofPointToTrackdZ[iSt].size(); i++)
+ TofPointToTrackdZ[iSt][i] = 100000;
- for (UInt_t iMC = 0; iMC < ToFPointsMatch.size(); iMC++) {
- if (ToFPointsMatch[iMC] == 0) continue;
+ for (Int_t iMC = 0; iMC < fTofPoints->Size(iFile, iEvent); iMC++) {
- int eventNr = ToFPointsMatch[iMC]->GetEntry();
+ CbmL1MCPoint MC;
- if (eventNr != iEvent) continue;
- CbmL1MCPoint MC;
+ if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 4)) {
- if (!ReadMCPoint(&MC, ToFPointsMatch[iMC]->GetIndex(), ToFPointsMatch[iMC]->GetFile(),
- ToFPointsMatch[iMC]->GetEntry(), 4)) {
+ Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
+ DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
+ if (trk_it == dFEI2vMCTracks.end()) continue;
+ Int_t IND_Track = trk_it->second;
MC.iStation = -1;
L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations + NTrdStations;
@@ -360,29 +388,46 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
MC.iStation = (MC.z > sta[iSt].z[0] - 6)
? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
: MC.iStation;
+ if (MC.iStation < 0) continue;
+ assert(MC.iStation >= 0);
+ int iTofSta = MC.iStation - (NMvdStations + NStsStations + NMuchStations + NTrdStations);
+ float dz = TofPointToTrackdZ[iTofSta][IND_Track];
+ if (MC.iStation >= 0)
+ if (fabs(sta[MC.iStation].z[0] - MC.z) < dz) TofPointToTrack[iTofSta][IND_Track] = iMC;
+ TofPointToTrackdZ[iTofSta][IND_Track] = fabs(sta[MC.iStation].z[0] - MC.z);
+ }
+ }
- Double_t dtrck = dFEI(iFile, eventNr, MC.ID);
- DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
- if (trk_it == dFEI2vMCTracks.end()) continue;
+ for (int iTofSta = 0; iTofSta < NTOFStation; iTofSta++)
+ for (unsigned int iMC = 0; iMC < TofPointToTrack[iTofSta].size(); iMC++) {
- Int_t IND_Track = trk_it->second;
+ if (TofPointToTrack[iTofSta][iMC] == 0) continue;
- vMCTracks[IND_Track].Points.push_back_no_warning(vMCPoints.size());
+ CbmL1MCPoint MC;
- MC.ID = trk_it->second;
+ if (!ReadMCPoint(&MC, TofPointToTrack[iTofSta][iMC], iFile, iEvent, 4)) {
- vMCPoints.push_back(MC);
- vMCPoints_in_Time_Slice.push_back(0);
+ MC.iStation = -1;
+ L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations + NTrdStations;
+ for (Int_t iSt = 0; iSt < NTOFStation; iSt++)
+ MC.iStation = (MC.z > sta[iSt].z[0] - 5)
+ ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
+ : MC.iStation;
+
+ vMCTracks[iMC].Points.push_back_no_warning(vMCPoints.size());
- dFEI2vMCPoints.insert(DFEI2I::value_type(
- dFEI(iFile, eventNr,
- ToFPointsMatch[iMC]->GetIndex() + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints),
- vMCPoints.size() - 1));
- nTofPoints++;
+ MC.ID = iMC;
+
+ vMCPoints.push_back(MC);
+ vMCPoints_in_Time_Slice.push_back(0);
+
+ dFEI2vMCPoints.insert(DFEI2I::value_type(
+ dFEI(iFile, iEvent, TofPointToTrack[iTofSta][iMC] + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints),
+ vMCPoints.size() - 1));
+ nTofPoints++;
+ }
}
- }
}
-
} //time_slice
for (unsigned int iTr = 0; iTr < vMCTracks.size(); iTr++) {
@@ -462,7 +507,27 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // if listMvdHits
if (fVerbose >= 10) cout << "ReadEvent: mvd hits are gotten." << endl;
- if (listStsHits) {
+ if ((2 == fStsUseMcHit)) { // SG!! create TRD hits from TRD points
+
+ int firstDetStrip = NStrips;
+ for (int ip = firstStsPoint; ip < firstStsPoint + nStsPoints; ip++) {
+ const CbmL1MCPoint& p = vMCPoints[ip];
+// int mcTrack = p.ID;
+// if (mcTrack < 0) continue;
+// const CbmL1MCTrack& t = vMCTracks[mcTrack];
+ //if (t.p < 1) continue;
+ // if (t.Points.size() > 4) continue;
+ // cout << "sts mc: station " << p.iStation - NMvdStations << " x " << p.x << " y " << p.y << " z " << p.z << " t "
+ // << p.time << " mc " << p.ID << " p " << p.p << endl;
+ TmpHit th;
+ int DetId = 1;
+ th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+ tmpHits.push_back(th);
+ nStsHits++;
+ }
+ }
+
+ if (listStsHits && (2 != fStsUseMcHit)) {
Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : listStsHits->GetEntriesFast());
@@ -592,7 +657,28 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // for j
} // if listStsHits
- if (fUseMUCH && fMuchPixelHits) {
+ if ((2 == fMuchUseMcHit) && fUseMUCH) { // SG!! create TRD hits from TRD points
+ int firstDetStrip = NStrips;
+
+ for (int ip = firstMuchPoint; ip < firstMuchPoint + nMuchPoints; ip++) {
+ const CbmL1MCPoint& p = vMCPoints[ip];
+
+// int mcTrack = p.ID;
+// if (mcTrack < 0) continue;
+// const CbmL1MCTrack& t = vMCTracks[mcTrack];
+ //if (t.p < 1) continue;
+ // if (t.Points.size() > 4) continue;
+
+ TmpHit th;
+ int DetId = 2;
+ th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+
+ tmpHits.push_back(th);
+ nMuchHits++;
+ }
+ }
+
+ if (fUseMUCH && fMuchPixelHits && (2 != fMuchUseMcHit)) {
Int_t nEnt = fMuchPixelHits->GetEntriesFast();
@@ -627,7 +713,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
th.x = mh->GetX();
- th.y = mh->GetY() - 0.5;
+ th.y = mh->GetY();
th.z = mh->GetZ();
th.dx = mh->GetDx();
@@ -663,32 +749,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
if (trk_it == dFEI2vMCPoints.end()) continue;
th.iMC = trk_it->second;
- // th.track = vMCPoints[th.iMC].ID;
- // th.qp = vMCPoints[iMC].p;
- // if(matchHitMatch -> GetNofLinks() == 0) continue;
- // Float_t bestWeight = 0.f;
- // Float_t totalWeight = 0.f;
- // int iMCPoint = -1;
- // CbmLink link;
- //
- // CbmMuchPoint* pt = (CbmMuchPoint*) fMuchPoints->Get(
- // matchHitMatch->GetLink(0).GetFile(),
- // matchHitMatch->GetLink(0).GetEntry(),
- // matchHitMatch->GetLink(0).GetIndex());
- // // double mom = sqrt(pt->GetPxOut()*pt->GetPxOut()+pt->GetPyOut()*pt->GetPyOut()+pt->GetPzOut()*pt->GetPzOut());
- // // th.p = mom;
- // // th.q = pt->GetTrackID();//(L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(link.GetFile(),link.GetEntry(), pt->GetTrackID()) ))->GetCharge();
- //
- // th.x = pt->GetX( th.z );// + gRandom->Gaus(0,th.dx);
- //
- // th.y = pt->GetY(th.z);// + gRandom->Gaus(0,th.dy);
- // th.time = pt->GetTime(); //+ gRandom->Gaus(0,th.dt);
- //
- // L1Station& st = algo->vStations[th.iStation];
- // th.u_front =
- // th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- // th.u_back =
- // th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ if ((1 == fMuchUseMcHit) && (th.iMC > -1))
+ th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
}
}
}
@@ -700,50 +762,26 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // for j
} // if listMvdHits
- if (0 && (fTrackingMode == L1Algo::TrackingMode::kGlobal) && fUseTRD) { // SG!! create TRD hits from TRD points
+ if ((2 == fTrdUseMcHit) && fTrdHitMatches && fUseTRD) { // SG!! create TRD hits from TRD points
int firstDetStrip = NStrips;
for (int ip = firstTrdPoint; ip < firstTrdPoint + nTrdPoints; ip++) {
const CbmL1MCPoint& p = vMCPoints[ip];
- int mcTrack = p.ID;
- if (mcTrack < 0) continue;
- const CbmL1MCTrack& t = vMCTracks[mcTrack];
- //if (t.p < 1) continue;
- if (t.Points.size() > 4) continue;
+// int mcTrack = p.ID;
+// if (mcTrack < 0) continue;
+// const CbmL1MCTrack& t = vMCTracks[mcTrack];
- cout << "trd mc: station " << p.iStation - NMvdStations - NStsStations - NMuchStations << " x " << p.x << " y "
- << p.y << " z " << p.z << " t " << p.time << " mc " << p.ID << " p " << p.p << endl;
TmpHit th;
- th.ExtIndex = 0;
- th.Det = 3;
- th.id = tmpHits.size();
- th.iStation = p.iStation;
- th.time = p.time;
- th.dt = 1.e4;
- th.iStripF = firstDetStrip + ip;
- th.iStripB = th.iStripF;
- if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
- th.x = p.x;
- th.y = p.y;
- th.z = p.z;
- float sigma = 1.;
- th.dx = sigma;
- th.dy = sigma;
- th.dxy = 0;
- th.du = sigma;
- th.dv = sigma;
- L1Station& st = algo->vStations[th.iStation];
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
- th.iMC = ip;
+ int DetId = 3;
+ th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
tmpHits.push_back(th);
nTrdHits++;
}
}
- if (1 && fUseTRD && listTrdHits) {
+ if (fUseTRD && listTrdHits && (2 != fTrdUseMcHit)) {
int firstDetStrip = NStrips;
vector<bool> mcUsed(nTrdPoints, 0);
@@ -799,7 +837,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.iMC = -1;
int iMC = -1;
- if (1 && fPerformance && fTrdHitMatches) {
+ if (fPerformance && fTrdHitMatches) {
CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j));
@@ -809,7 +847,11 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
// th.track = vMCPoints[th.iMC].ID;
- if (fTrackingMode == L1Algo::TrackingMode::kGlobal) { //SG!!! replace hits by MC points
+ if ((1 == fTrdUseMcHit) && (th.iMC > -1)) //SG!!! replace hits by MC points
+
+ th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
+
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { //SG!!! replace hits by MC points
CbmTrdPoint* pt =
(CbmTrdPoint*) fTrdPoints->Get(trdHitMatch->GetLink(0).GetFile(), trdHitMatch->GetLink(0).GetEntry(),
@@ -854,8 +896,30 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
} // for j
} // if listTrdHits
+ if ((2 == fTofUseMcHit) && fUseTOF) { // SG!! create TRD hits from TRD points
+ int firstDetStrip = NStrips;
+
+ for (int ip = firstTofPoint; ip < firstTofPoint + nTofPoints; ip++) {
+ const CbmL1MCPoint& p = vMCPoints[ip];
+
+// int mcTrack = p.ID;
+ // if (mcTrack < 0) continue;
+ //const CbmL1MCTrack& t = vMCTracks[mcTrack];
+ //if (t.p < 1) continue;
+ //if (t.Points.size() > 4) continue;
+
+ TmpHit th;
+
+ int DetId = 4;
+ th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+ tmpHits.push_back(th);
+
+ nTofHits++;
+ }
+ }
- if (fUseTOF && fTofHits) {
+
+ if (fUseTOF && fTofHits && (2 != fTofUseMcHit)) {
int firstDetStrip = NStrips;
@@ -924,36 +988,25 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (fPerformance) {
- if (ToFPointsMatch[j] == 0) continue;
-
- // th.iMC = j+nMvdPoints+nStsPoints+nTrdPoints+nMuchPoints;
-
- Int_t iFile = ToFPointsMatch[j]->GetFile();
- Int_t iEvent = ToFPointsMatch[j]->GetEntry();
-
- Int_t iIndex = ToFPointsMatch[j]->GetIndex() + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints;
+ CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(j));
- Double_t dtrck = dFEI(iFile, iEvent, iIndex);
- DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
- if (trk_it == dFEI2vMCPoints.end()) continue;
- th.iMC = trk_it->second;
- // th.track = vMCPoints[th.iMC].ID;
- // th.qp = vMCPoints[th.iMC].p;
- // CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(ToFPointsMatch[j]->GetFile(),ToFPointsMatch[j]->GetEntry(),ToFPointsMatch[j]->GetIndex());
+ if (0 < matchHitMatch->GetNofLinks()) {
+ Int_t iMC = matchHitMatch->GetLink(0).GetIndex();
+ Int_t iFile = matchHitMatch->GetLink(0).GetFile();
+ Int_t iEvent = matchHitMatch->GetLink(0).GetEntry();
+ CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fTofPoints->Get(iFile, iEvent, iMC));
+ pt->GetTrackID();
+ Double_t dtrck = dFEI(iFile, iEvent, pt->GetTrackID());
+ DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
+ if (trk_it == dFEI2vMCPoints.end()) continue;
+ th.iMC = TofPointToTrack[sttof][trk_it->second];
- // th.x = pt->GetX() + gRandom->Gaus(0,th.dx);
- //
- // th.y = pt->GetY()+ gRandom->Gaus(0,th.dy);
- // th.time = pt->GetTime()+ gRandom->Gaus(0,th.dt);
- //
- // L1Station &st = algo->vStations[th.iStation];
- // th.u_front = th.x*st.frontInfo.cos_phi[0] + th.y*st.frontInfo.sin_phi[0];
- // th.u_back = th.x* st.backInfo.cos_phi[0] + th.y*st.backInfo.sin_phi[0];
+ if ((1 == fTofUseMcHit) && (th.iMC > -1)) th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
+ }
}
-
tmpHits.push_back(th);
nTofHits++;
@@ -1200,9 +1253,11 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int M
Double_t StartTime = fTimeSlice->GetStartTime();
Double_t EndTime = fTimeSlice->GetEndTime();
Double_t Time_MC_point = pt->GetTime() + fEventList->GetEventTime(event, file);
- if (Time_MC_point < StartTime) return 1;
-
- if (Time_MC_point > EndTime) return 1;
+
+ cout<<StartTime<<" StartTime "<<EndTime<<" EndTime "<<Time_MC_point<<" Time_MC_point "<<fEventList->GetEventTime(event, file)<<endl;
+ // if (Time_MC_point < StartTime) return 1;
+ //
+ // if (Time_MC_point > EndTime) return 1;
}
pt->Position(xyzI);
@@ -1242,12 +1297,12 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int M
CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fTofPoints->Get(file, event, iPoint)); // file, event, object
if (!pt) return 1;
if (!fLegacyEventMode) {
- Double_t StartTime = fTimeSlice->GetStartTime();
- Double_t EndTime = fTimeSlice->GetEndTime();
- Double_t Time_MC_point = pt->GetTime() + fEventList->GetEventTime(event, file);
- if (Time_MC_point < StartTime) return 1;
-
- if (Time_MC_point > EndTime) return 1;
+ // Double_t StartTime = fTimeSlice->GetStartTime();
+ // Double_t EndTime = fTimeSlice->GetEndTime();
+ // Double_t Time_MC_point = pt->GetTime() + fEventList->GetEventTime(event, file);
+ // if (Time_MC_point < StartTime) return 1;
+ //
+ // if (Time_MC_point > EndTime) return 1; continue!!!!
}
pt->Position(xyzI);
@@ -1312,7 +1367,7 @@ void CbmL1::HitMatch()
for (int iH = 0; iH < NHits; iH++) {
CbmL1Hit& hit = vStsHits[iH];
- if (hit.Det == 1) {
+ if ((hit.Det == 1) && (2 != fStsUseMcHit)) {
CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(vStsHits[iH].extIndex));
int iP = -1;
@@ -1395,7 +1450,7 @@ void CbmL1::HitMatch()
} // if no clusters
}
- if (hit.Det != 1) { // the hit is not from STS
+ if ((hit.Det != 1) || (2 == fStsUseMcHit)) { // the hit is not from STS
int iP = vHitMCRef[iH];
if (iP >= 0) {
hit.mcPointIds.push_back_no_warning(iP);
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 97f9be3d83..5b2a9dbfb0 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -132,6 +132,7 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
double f_sigma = geo[ind++];
double b_phi = geo[ind++];
double b_sigma = geo[ind++];
+ double dt = geo[ind++];
double c_f = cos(f_phi);
double s_f = sin(f_phi);
double c_b = cos(b_phi);
@@ -144,6 +145,7 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
st.backInfo.cos_phi = c_b;
st.backInfo.sin_phi = s_b;
st.backInfo.sigma2 = b_sigma * b_sigma;
+ st.dt = dt;
if (fabs(b_phi - f_phi) < .1) {
double th = b_phi - f_phi;
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index db29a74cac..eb8acd6b82 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -12,7 +12,6 @@
class L1Station {
public:
-
// TODO: test this constructors
//L1Station(const L1Station &) = default;
//L1Station & operator=(const L1Station &) = default;
@@ -25,6 +24,7 @@ public:
fvec Rmin {0};
fvec Rmax {0};
fvec Sy {0};
+ fvec dt {0};
L1MaterialInfo materialInfo {};
L1FieldSlice fieldSlice {};
L1UMeasurementInfo frontInfo {};
--
GitLab