CbmRoot
CbmTrdToTofVector.cxx
Go to the documentation of this file.
1 
5 #include "CbmTrdToTofVector.h"
6 #include "CbmDefs.h"
7 #include "CbmMCDataArray.h"
8 #include "CbmMCDataManager.h"
9 #include "CbmMatch.h"
10 #include "CbmMuchTrack.h"
11 #include "CbmSetup.h"
12 #include "CbmTofAddress.h"
13 #include "CbmTofDigi.h"
14 #include "CbmTofHit.h"
15 #include "CbmTofPoint.h"
16 
17 #include "FairRootManager.h"
18 
19 #include <TClonesArray.h>
20 #include <TGeoBBox.h>
21 #include <TGeoManager.h>
22 #include <TMath.h>
23 #include <TMatrixDSym.h>
24 #include <TVectorD.h>
25 
26 #include <iostream>
27 
28 using std::cout;
29 using std::endl;
30 using std::map;
31 using std::multimap;
32 using std::pair;
33 using std::set;
34 
35 // ----- Default constructor -------------------------------------------
36 CbmTrdToTofVector::CbmTrdToTofVector()
37  : FairTask("Trd2TofVec")
38  , fTrackArray(NULL)
39  , fNofTracks(0)
40  , fHits(NULL)
41  , fHitMatches(NULL)
42  , fPoints(NULL)
43  , fDigis(NULL)
44  , fDigiMatches(NULL)
45  , fTrdTracks(NULL)
46  , fErrX(1.0)
47  , fErrY(1.0)
48  , fCutChi2(24.0) // chi2/ndf = 6 for 4 hits
49 {}
50 // -------------------------------------------------------------------------
51 
52 // ----- Destructor ----------------------------------------------------
53 CbmTrdToTofVector::~CbmTrdToTofVector() {}
54 // -------------------------------------------------------------------------
55 
56 // ----- Public method Init (abstract in base class) --------------------
57 InitStatus CbmTrdToTofVector::Init() {
58 
59  // Get and check FairRootManager
60  FairRootManager* ioman = FairRootManager::Instance();
61  if (ioman == NULL) Fatal("Init", "RootManager not instantiated!");
62 
63  // Create and register TrdTrack array (if necessary)
64  //fTrackArray = static_cast<TClonesArray*> (ioman->GetObject("MuchVector"));
65  //if (fTrackArray == NULL) {
66  fTrackArray = new TClonesArray("CbmMuchTrack", 100);
67  ioman->Register("TofVector", "Tof", fTrackArray, kTRUE);
68  //}
69  fHits = static_cast<TClonesArray*>(ioman->GetObject("TofHit"));
70  fHitMatches = static_cast<TClonesArray*>(ioman->GetObject("TofDigiMatch"));
71  //fPoints = static_cast<TClonesArray*> (ioman->GetObject("TofPoint"));
72  CbmMCDataManager* mcManager =
73  (CbmMCDataManager*) ioman->GetObject("MCDataManager");
74  if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
75  fPoints = mcManager->InitBranch("TofPoint");
76  fDigis = static_cast<TClonesArray*>(ioman->GetObject("TofDigi"));
77  fDigiMatches =
78  static_cast<TClonesArray*>(ioman->GetObject("TofDigiMatchPoints"));
79  fTrdTracks = static_cast<TClonesArray*>(ioman->GetObject("TrdVector"));
80 
81  TString tag;
82  CbmSetup::Instance()->GetGeoTag(kTof, tag);
83  cout << " ******* TOF tag ******** " << tag << endl;
84  tag.Prepend("tof_");
85 
86  //CbmTofGeoHandler geoHandler;
87  //geoHandler.Init();
88 
89  // Get TOF location in Z
90  TGeoVolume* tofV = gGeoManager->GetVolume(tag);
91  TGeoBBox* shape = (TGeoBBox*) tofV->GetShape();
92  shape->GetAxisRange(3, fZ[0], fZ[1]);
93  cout << " TOF span in Z: " << fZ[0] << " " << fZ[1] << endl;
94 
95  return kSUCCESS;
96 }
97 // -------------------------------------------------------------------------
98 
99 // ----- SetParContainers -------------------------------------------------
100 void CbmTrdToTofVector::SetParContainers() {
101  /*
102  fDigiPar = (CbmTrdDigiPar*) FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdDigiPar");
103  cout << " ******* digiPar ******** " << fDigiPar << endl;
104  exit(0);
105  */
106 }
107 // -------------------------------------------------------------------------
108 
109 // ----- Public method Exec --------------------------------------------
110 void CbmTrdToTofVector::Exec(Option_t* opt) {
111 
112  fTrackArray->Delete();
113 
114  // Do all processing
115 
116  // Get hits
117  GetHits();
118 
119  // Build vectors (extrapolate TRD vectors thru the thickness of TOF)
120  MakeVectors();
121 
122  // Remove vectors with wrong orientation
123  // (using empirical cuts for omega muons at 8 GeV)
124  //CheckParams();
125 
126  // Match TOF hits to TRD vectors
127  MatchTofToTrd();
128 
129  // Remove clones
130  //RemoveClones();
131 
132  // Store vectors
133  StoreVectors();
134 }
135 // -------------------------------------------------------------------------
136 
137 // ----- Public method Finish ------------------------------------------
138 void CbmTrdToTofVector::Finish() {}
139 // -------------------------------------------------------------------------
140 
141 // ----- Private method GetHits ----------------------------------------
142 void CbmTrdToTofVector::GetHits() {
143  // Arrange hits according to their Phi-angle
144 
145  fHitX.clear();
146  fHitY.clear();
147  fHitIds.clear();
148  fHitTime.clear();
149 
150  Int_t nHits = fHits->GetEntriesFast(), sel = 0;
151 
152  for (Int_t i = 0; i < nHits; ++i) {
153  CbmTofHit* hit = (CbmTofHit*) fHits->UncheckedAt(i);
154 
156  sel = 1; //SelectHitId(hit);
157  if (!sel) continue;
158  //
159  fHitX.insert(pair<Double_t, Int_t>(hit->GetX(), i));
160  fHitY.insert(pair<Double_t, Int_t>(hit->GetY(), i));
161  CbmMatch* match = (CbmMatch*) fHitMatches->UncheckedAt(i);
162  Int_t nlinks = match->GetNofLinks();
163  set<Int_t> ids;
164  fHitIds[i] = ids;
165 
166  for (Int_t il = 0; il < nlinks; ++il) {
167  const CbmLink link = match->GetLink(il);
168  CbmTofDigi* digi = (CbmTofDigi*) fDigis->UncheckedAt(link.GetIndex());
169  //AZ CbmMatch *digiM = (CbmMatch*) fDigiMatches->UncheckedAt(link.GetIndex());
170  CbmMatch* digiM = digi->GetMatch();
171  Int_t npoints = digiM->GetNofLinks();
172 
173  for (Int_t ip = 0; ip < npoints; ++ip) {
174  const CbmLink link1 = digiM->GetLink(ip);
175  //CbmTofPoint *point = (CbmTofPoint*) fPoints->UncheckedAt(link1.GetIndex());
176  CbmTofPoint* point = (CbmTofPoint*) fPoints->Get(
177  link1.GetFile(), link1.GetEntry(), link1.GetIndex());
178  fHitIds[i].insert(point->GetTrackID());
179  if (fHitTime.find(i) == fHitTime.end())
180  fHitTime[i] = point->GetTime();
181  else
182  fHitTime[i] = TMath::Min(point->GetTime(), fHitTime[i]);
183  //cout << point->GetTrackID() << endl;
184  }
185  }
186  }
187 }
188 // -------------------------------------------------------------------------
189 
190 // ----- Private method MakeVectors ------------------------------------
191 void CbmTrdToTofVector::MakeVectors() {
192  // Make vectors (extrapolate TRD vectors through the TOF thickness)
193 
194  Int_t nvec = fVectors.size();
195  for (Int_t j = 0; j < nvec; ++j)
196  delete fVectors[j];
197  fVectors.clear();
198 
199  fLineX.clear();
200 
201  Int_t nTrd = fTrdTracks->GetEntriesFast();
202  cout << " TRD tracks: " << nTrd << endl;
203 
204  for (Int_t iv = 0; iv < nTrd; ++iv) {
205  CbmMuchTrack* tr = (CbmMuchTrack*) fTrdTracks->UncheckedAt(iv);
206  const FairTrackParam& param = *tr->GetParamLast();
207  TLine line;
208  line.SetUniqueID(iv);
209 
210  for (Int_t i = 0; i < 2; ++i) {
211  Double_t dz = fZ[i] - param.GetZ();
212  if (i == 0) {
213  line.SetX1(param.GetX() + param.GetTx() * dz);
214  line.SetY1(param.GetY() + param.GetTy() * dz);
215  } else {
216  line.SetX2(param.GetX() + param.GetTx() * dz);
217  line.SetY2(param.GetY() + param.GetTy() * dz);
218  }
219  }
220  fLineX.insert(
221  pair<Double_t, TLine>(TMath::Min(line.GetX1(), line.GetX2()), line));
222  }
223 }
224 // -------------------------------------------------------------------------
225 
226 // ----- Private method MatchTofToTrd ----------------------------------
227 void CbmTrdToTofVector::MatchTofToTrd() {
228  // Match TOF hits to TRD vectors
229 
230  //const Double_t window = 30.0;
231  const Double_t window = 50.0;
232  multimap<Double_t, Int_t>::iterator mitb, mite, mit1;
233  multimap<Double_t, pair<Int_t, Int_t>> rads;
234  map<pair<Int_t, Int_t>, pair<Double_t, Double_t>> dtdl;
235  set<Int_t> setVec, setHit;
236 
237  for (multimap<Double_t, TLine>::iterator mit = fLineX.begin();
238  mit != fLineX.end();
239  ++mit) {
240  Double_t rmin = 999999, dtmin = 0, dlmin = 0;
241  TVector3 trd(mit->second.GetX2() - mit->second.GetX1(),
242  mit->second.GetY2() - mit->second.GetY1(),
243  0.0);
244  Double_t trdLeng = trd.Mag();
245  Int_t indVec = mit->second.GetUniqueID();
246  CbmMuchTrack* tr = (CbmMuchTrack*) fTrdTracks->UncheckedAt(indVec);
247  Int_t id = tr->GetFlag(), ihit = -1;
248 
249  // Select TOF hits for matching (apply coordinate windows)
250  mitb = fHitX.lower_bound(mit->first - window);
251  mite = fHitX.upper_bound(
252  TMath::Max(mit->second.GetX1(), mit->second.GetX2()) + window);
253  set<Int_t> inds;
254  for (mit1 = mitb; mit1 != mite; ++mit1)
255  inds.insert(mit1->second);
256  mitb = fHitY.lower_bound(
257  TMath::Min(mit->second.GetY1(), mit->second.GetY2()) - window);
258  mite = fHitY.upper_bound(
259  TMath::Max(mit->second.GetY1(), mit->second.GetY2()) + window);
260 
261  for (mit1 = mitb; mit1 != mite; ++mit1) {
262  if (inds.find(mit1->second) == inds.end()) continue; // outside window
263  CbmTofHit* hit = (CbmTofHit*) fHits->UncheckedAt(mit1->second);
264  TVector3 tof(hit->GetX() - mit->second.GetX1(),
265  hit->GetY() - mit->second.GetY1(),
266  0.0);
267  Double_t dt = TMath::Abs(trd.Cross(tof).Z()) / trdLeng;
268  Double_t dl = trd * tof / trdLeng;
269  if (dl > trdLeng)
270  dl -= trdLeng;
271  else if (dl > 0)
272  dl = 0;
273  Double_t rad = dl * dl + dt * dt;
274  // Save matches
275  rads.insert(pair<Double_t, pair<Int_t, Int_t>>(
276  rad, pair<Int_t, Int_t>(indVec, mit1->second)));
277  dtdl[pair<Int_t, Int_t>(indVec, mit1->second)] =
278  pair<Double_t, Double_t>(dt, dl);
279  setVec.insert(indVec);
280  setHit.insert(mit1->second);
281  /*
282  if (rad < rmin) {
283  rmin = rad;
284  dtmin = dt;
285  dlmin = dl;
286  ihit = mit1->second;
287  }
288  */
289  }
290  /*
291  if (ihit >= 0) {
292  Bool_t match = kFALSE;
293  if (fHitIds[ihit].find(id) != fHitIds[ihit].end()) match = kTRUE;
294  //cout << dtmin << " " << dlmin << " " << id << " " << match << endl;
295  CbmMuchTrack *trTof = new CbmMuchTrack(*tr);
296  trTof->SetPreviousTrackId(mit->second.GetUniqueID()); // index of TRD vector
297  FairTrackParam param;
298  param.SetX(dtmin); // use storage
299  param.SetY(dlmin);
300  trTof->SetParamLast(&param);
301  trTof->SetFlag(match);
302  trTof->SetUniqueID(id);
303  fVectors.push_back(trTof);
304  }
305  */
306  } // for (multimap<Double_t,TLine>::iterator mit = fLineX.begin();
307 
308  // Create vectors
309  for (multimap<Double_t, pair<Int_t, Int_t>>::iterator rmit = rads.begin();
310  rmit != rads.end();
311  ++rmit) {
312  Int_t indVec = rmit->second.first, ihit = rmit->second.second;
313  if (setVec.find(indVec) == setVec.end())
314  continue; // already matched vector
315  if (setHit.find(ihit) == setHit.end()) continue; // already matched hit
316 
317  Bool_t match = kFALSE;
318  CbmMuchTrack* tr = (CbmMuchTrack*) fTrdTracks->UncheckedAt(indVec);
319  Int_t id = tr->GetFlag();
320  if (fHitIds[ihit].find(id) != fHitIds[ihit].end()) match = kTRUE;
321  //cout << dtmin << " " << dlmin << " " << id << " " << match << endl;
322  CbmMuchTrack* trTof = new CbmMuchTrack(*tr);
323  trTof->SetPreviousTrackId(indVec); // index of TRD vector
324  FairTrackParam param;
325  param.SetX(dtdl[rmit->second].first); // use storage
326  param.SetY(dtdl[rmit->second].second);
327  param.SetZ(rmit->first);
328  trTof->SetParamLast(&param);
329  trTof->SetFlag(match);
330  trTof->SetUniqueID(id);
331  trTof->SetNDF(ihit); // TOF hit index - use storage
332  trTof->SetChiSq(fHitTime[ihit]); // min. time - use storage
333  fVectors.push_back(trTof);
334 
335  setVec.erase(indVec);
336  if (setVec.size() == 0) break;
337  setHit.erase(ihit);
338  if (setHit.size() == 0) break;
339  }
340 }
341 // -------------------------------------------------------------------------
342 
343 // ----- Private method CheckParams ------------------------------------
344 void CbmTrdToTofVector::CheckParams() {
345  // Remove vectors with wrong orientation
346  // using empirical cuts for omega muons at 8 Gev
347 
348  const Double_t cut[2] = {0.6, 0.6}; // !!! empirical !!!
349 
350  Int_t nvec = fVectors.size();
351 
352  for (Int_t iv = 0; iv < nvec; ++iv) {
353  CbmMuchTrack* vec = fVectors[iv];
354  const FairTrackParam* params = vec->GetParamFirst();
355  Double_t dTx = params->GetTx() - params->GetX() / params->GetZ();
356  if (TMath::Abs(dTx) > cut[0])
357  vec->SetChiSq(-1.0);
358  else {
359  Double_t dTy = params->GetTy() - params->GetY() / params->GetZ();
360  if (TMath::Abs(dTy) > cut[1]) vec->SetChiSq(-1.0);
361  }
362  }
363 
364  for (Int_t iv = nvec - 1; iv >= 0; --iv) {
365  CbmMuchTrack* vec = fVectors[iv];
366  if (vec->GetChiSq() < 0) {
367  delete fVectors[iv];
368  fVectors.erase(fVectors.begin() + iv);
369  }
370  }
371  cout << " Vectors after parameter check: " << nvec << " " << fVectors.size()
372  << endl;
373 }
374 // -------------------------------------------------------------------------
375 
376 // ----- Private method RemoveClones -----------------------------------
377 void CbmTrdToTofVector::RemoveClones() {
378  // Remove clone vectors (having at least 3 the same hits)
379 
380  //Int_t nthr = 3, planes[20];
381  Int_t nthr = 2, planes[20];
382 
383  Int_t nvec = fVectors.size();
384 
385  // Do sorting according to "quality"
386  multimap<Double_t, CbmMuchTrack*> qMap;
387  multimap<Double_t, CbmMuchTrack*>::iterator it, it1;
388 
389  for (Int_t i = 0; i < nvec; ++i) {
390  CbmMuchTrack* vec = fVectors[i];
391  Double_t qual =
392  vec->GetNofHits() + (99 - TMath::Min(vec->GetChiSq(), 99.0)) / 100;
393  qMap.insert(pair<Double_t, CbmMuchTrack*>(-qual, vec));
394  }
395 
396  for (it = qMap.begin(); it != qMap.end(); ++it) {
397  CbmMuchTrack* vec = it->second;
398  if (vec->GetChiSq() < 0) continue;
399  for (Int_t j = 0; j < fgkPlanes; ++j)
400  planes[j] = -1;
401 
402  Int_t nhits = vec->GetNofHits();
403  for (Int_t ih = 0; ih < nhits; ++ih) {
404  CbmTofHit* hit = (CbmTofHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
405  Int_t lay = hit->GetPlaneId();
406  planes[lay] = vec->GetHitIndex(ih);
407  }
408 
409  it1 = it;
410  for (++it1; it1 != qMap.end(); ++it1) {
411  CbmMuchTrack* vec1 = it1->second;
412  if (vec1->GetChiSq() < 0) continue;
413  Int_t nsame = 0, same[fgkPlanes] = {0};
414 
415  Int_t nhits1 = vec1->GetNofHits();
416  //nthr = TMath::Min(nhits,nhits1) / 2;
417  //nthr = TMath::Min(nhits,nhits1) * 0.75;
418  for (Int_t ih = 0; ih < nhits1; ++ih) {
419  CbmTofHit* hit = (CbmTofHit*) fHits->UncheckedAt(vec1->GetHitIndex(ih));
420  Int_t lay = hit->GetPlaneId();
421  if (planes[lay] >= 0) {
422  if (vec1->GetHitIndex(ih) == planes[lay]) same[lay] = 1;
423  //else same[lay] = 0;
424  }
425  }
426  for (Int_t lay = 0; lay < fgkPlanes; ++lay)
427  nsame += same[lay];
428 
429  if (nsame >= nthr) {
430  // Too many the same hits
431  Int_t clone = 0;
432  if (nhits > nhits1 + 0)
433  clone = 1;
434  else if (vec->GetChiSq() * 1 <= vec1->GetChiSq())
435  clone = 1; // the same number of hits on 2 tracks
436  //else if (vec->GetChiSq() * 1.5 <= vec1->GetChiSq()) clone = 1; // the same number of hits on 2 tracks
437  if (clone) vec1->SetChiSq(-1.0);
438  }
439  }
440  } // for (it = qMap.begin();
441 
442  for (Int_t iv = nvec - 1; iv >= 0; --iv) {
443  CbmMuchTrack* vec = fVectors[iv];
444  if (vec->GetChiSq() < 0) {
445  delete fVectors[iv];
446  fVectors.erase(fVectors.begin() + iv);
447  }
448  }
449  cout << " CbmTrdToTofVector:: Vectors after clones removed: " << nvec << " "
450  << fVectors.size() << endl;
451 }
452 // -------------------------------------------------------------------------
453 
454 // ----- Private method StoreVectors -----------------------------------
455 void CbmTrdToTofVector::StoreVectors() {
456  // Store vectors (CbmMuchTracks) into TClonesArray
457 
458  Int_t ntrs = fTrackArray->GetEntriesFast();
459  Int_t nHitsTot = fHits->GetEntriesFast();
460 
461  //set<Int_t> usedHits;
462  Int_t nvec = fVectors.size();
463 
464  for (Int_t iv = 0; iv < nvec; ++iv) {
465  CbmMuchTrack* tr =
466  new ((*fTrackArray)[ntrs++]) CbmMuchTrack(*(fVectors[iv]));
467  //cout << " Track: " << tr->GetNofHits() << endl;
468  //for (Int_t j = 0; j < tr->GetNofHits(); ++j) cout << j << " " << tr->GetHitIndex(j) << " " << fVectors[ist][iv]->GetHitIndex(j) << endl;
469  // Set hit flag (to check Lit tracking)
470  //Int_t nhits = tr->GetNofHits();
471 
472  // Save IDs of contributing points to the hit
473  CbmMatch* match = new CbmMatch;
474  Int_t ihit = fVectors[iv]->GetNDF();
475  for (set<Int_t>::iterator sit = fHitIds[ihit].begin();
476  sit != fHitIds[ihit].end();
477  ++sit)
478  match->AddLink(1.0, *sit);
479  tr->SetMatch(match);
480  }
481 }
482 // -------------------------------------------------------------------------
483 
484 // ----- Private method ProcessPlane -----------------------------------
485 void CbmTrdToTofVector::ProcessPlane(Int_t lay, Int_t patt, Int_t flag0) {
486  // Main processing engine (recursively adds layer hits to the vector)
487 
488  //const Double_t cut[2] = {0.8, 0.8}; // !!! empirical !!!
489  const Double_t cut[2] = {0.6, 0.6}; // !!! empirical !!!
490 
491  Double_t pars[4] = {0.0};
492  Int_t flag = 0;
493  multimap<Int_t, Int_t>::iterator mit;
494  // Clear bits
495  patt &= ~(1 << lay);
496  /*
497  for (mit = fHitPl[lay].begin(); mit != fHitPl[lay].end(); ++mit) {
498  Int_t indx = mit->second;
499 
500  CbmTofHit *hit = (CbmTofHit*) fHits->UncheckedAt(indx);
501  fXy[lay][0] = hit->GetX();
502  fXy[lay][1] = hit->GetY();
503  fXy[lay][2] = hit->GetDx();
504  fXy[lay][3] = hit->GetDy();
505  //fXy[lay][2] = TMath::Max (hit->GetDx(),fErrX[ista]);
506  //fXy[lay][3] = TMath::Max (hit->GetDy(),fErrY[ista]);
507  fXy[lay][4] = hit->GetZ();
508  // Clear bits
509  patt &= ~(1 << lay);
510 
511  // Check slopes
512  Int_t lay0 = 0;
513  if (!(patt & (1 << lay0))) ++lay0;;
514  Double_t dx = fXy[lay][0] - fXy[lay0][0];
515  Double_t dz = fXy[lay][4] - fXy[lay0][4];
516  Double_t dTx = TMath::Abs(dx/dz) - TMath::Abs(fXy[lay][0]/fXy[lay][4]);
517  if (TMath::Abs(dTx) > cut[0]) continue;
518  Double_t dy = fXy[lay][1] - fXy[lay0][1];
519  Double_t dTy = TMath::Abs(dy/dz) - TMath::Abs(fXy[lay][1]/fXy[lay][4]);
520  if (TMath::Abs(dTy) > cut[1]) continue;
521 
522  fXyi[lay][1] = indx;
523  flag = 1;
524 
525  // Set bit
526  patt |= (1 << lay);
527  //cout << lay << " " << patt << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << " " << endl;
528 
529  if (lay < fgkPlanes - 1) ProcessPlane(lay+1, patt, flag);
530  else if (TString::Itoa(patt,2).CountChar('1') > 2) {
531  // Straight line fit of the vector with > 2 hits
532  FindLine(patt, pars);
533  Double_t chi2 = FindChi2(patt, pars);
534  cout << " *** Stat: " << chi2 << " " << pars[0] << " " << pars[1] << endl;
535  if (chi2 <= fCutChi2) AddVector(patt, chi2, pars); // add vector to the temporary container
536  }
537  } // for (mit = fHitPl[lay].begin();
538 
539  //if (flag == 0) {
540  if (flag >= 0) {
541  // Missing plane
542  if (flag0 == 0) return; // 2 consecutive missing planes
543  // Clear bits
544  patt &= ~(1 << lay);
545  if (lay < fgkPlanes - 1) ProcessPlane(lay+1, patt, flag);
546  else if (TString::Itoa(patt,2).CountChar('1') > 2) {
547  // Straight line fit of the vector with > 2 hits
548  FindLine(patt, pars);
549  Double_t chi2 = FindChi2(patt, pars);
550  cout << " *** Stat: " << chi2 << " " << pars[0] << " " << pars[1] << endl;
551  if (chi2 <= fCutChi2) AddVector(patt, chi2, pars); // add vector to the temporary container
552  }
553  }
554  */
555 }
556 // -------------------------------------------------------------------------
557 
558 // ----- Private method FindLine ---------------------------------------
559 void CbmTrdToTofVector::FindLine(Int_t patt, Double_t* pars) {
560  // Fit of hits to the straight line
561 
562  // Solve system of linear equations
563  Bool_t ok = kFALSE, onoff;
564  TVectorD b(4);
565  /*
566  for (Int_t i = 0; i < fgkPlanes; ++i) {
567  onoff = patt & (1 << i);
568  if (!onoff) continue;
569  b[0] += fXy[i][0];
570  b[1] += fXy[i][1];
571  b[2] += fXy[i][0] * fDz[i];
572  b[3] += fXy[i][1] * fDz[i];
573  }
574 
575  //lu.Print();
576  TVectorD solve = fLus[patt]->Solve(b, ok);
577  //TVectorD solve = lu.TransSolve(b, ok);
578  //lu.Print();
579  for (Int_t i = 0; i < 4; ++i) pars[i] = solve[i];
580  //for (Int_t i = 0; i < 4; ++i) { cout << pars[i] << " "; if (i == 3) cout << endl; }
581  //Double_t y1 = cosa / sina * (uz[1][0] * cosa - uz[0][0]) + uz[1][0] * sina;
582  //Double_t y2 = -cosa / sina * (uz[2][0] * cosa - uz[0][0]) - uz[2][0] * sina;
583  //cout << " OK " << ok << " " << y1 << " " << y2 << endl;
584  */
585 }
586 // -------------------------------------------------------------------------
587 
588 // ----- Private method FindChi2 ---------------------------------------
589 Double_t CbmTrdToTofVector::FindChi2(Int_t patt, Double_t* pars) {
590  // Compute chi2 of the fit
591 
592  Double_t chi2 = 0, x = 0, y = 0;
593  Bool_t onoff;
594  /*
595  for (Int_t i = 0; i < fgkPlanes; ++i) {
596  onoff = patt & (1 << i);
597  if (!onoff) continue;
598  x = pars[0] + pars[2] * fDz[i];
599  y = pars[1] + pars[3] * fDz[i];
600  //Double_t dx = (x - fXy[i][0]) / fErrX;
601  //Double_t dy = (y - fXy[i][1]) / fErrY;
602  Double_t dx = (x - fXy[i][0]) / TMath::Max(fXy[i][2],fErrX);
603  Double_t dy = (y - fXy[i][1]) / TMath::Max(fXy[i][3],fErrY);
604  chi2 += dx * dx;
605  chi2 += dy * dy;
606  }
607  */
608  return chi2;
609 }
610 // -------------------------------------------------------------------------
611 
612 // ----- Private method AddVector --------------------------------------
613 void CbmTrdToTofVector::AddVector(Int_t patt, Double_t chi2, Double_t* pars) {
614  // Add vector to the temporary container (as a MuchTrack)
615 
616  Bool_t refit = kFALSE; //kTRUE;
617  TMatrixDSym cov(4);
618  /*
619  if (refit) {
620  // Refit line with individual hit errors
621  //cout << " Before: " << chi2 << endl;
622  //Refit(patt, chi2, pars, cov);
623  //cout << " After: " << chi2 << endl;
624  } else {
625  cov = *fMatr[patt];
626  cov *= (fErrX * fErrX); // the same error in X and Y
627  }
628  cov.ResizeTo(5,5);
629  cov(4,4) = 1.0;
630 
631  CbmMuchTrack *track = new CbmMuchTrack();
632  track->SetChiSq(chi2);
633  FairTrackParam par(pars[0], pars[1], fZ0, pars[2], pars[3], 0.0, cov);
634  track->SetParamFirst(&par);
635  par.SetZ(fZ0 + fDz[fgkPlanes-1]);
636  par.SetX(pars[0] + fDz[fgkPlanes-1] * pars[2]);
637  par.SetY(pars[1] + fDz[fgkPlanes-1] * pars[3]);
638  track->SetParamLast(&par);
639  //track->SetUniqueID(ista); // just to pass the value
640 
641  for (Int_t ipl = 0; ipl < fgkPlanes; ++ipl) {
642  if (!(patt & (1 << ipl))) continue;
643  track->AddHit(fXyi[ipl][1], kTRDHIT);
644  }
645  Int_t ndf = (track->GetNofHits() > 2) ? track->GetNofHits() * 2 - 4 : 1;
646  track->SetNDF(ndf);
647  SetTrackId(track); // set track ID as its flag
648  fVectors.push_back(track);
649  */
650 }
651 // -------------------------------------------------------------------------
652 
653 // ----- Private method SetTrackId -------------------------------------
654 void CbmTrdToTofVector::SetTrackId(CbmMuchTrack* vec) {
655  // Set vector ID as its flag (maximum track ID of its poins)
656 
657  map<Int_t, Int_t> ids;
658  Int_t nhits = vec->GetNofHits(), id = 0;
659 
660  /*
661  for (Int_t ih = 0; ih < nhits; ++ih) {
662  CbmTofHit *hit = (CbmTofHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
663  CbmTrdCluster *clus = (CbmTrdCluster*) fClusters->UncheckedAt(hit->GetRefId());
664  Int_t nDigis = clus->GetNofDigis();
665 
666  for (Int_t j = 0; j < nDigis; ++j) {
667  CbmMatch* digiM = (CbmMatch*) fDigiMatches->UncheckedAt(clus->GetDigi(j));
668  Int_t np = digiM->GetNofLinks();
669 
670  for (Int_t ip = 0; ip < np; ++ip) {
671  //CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(0).GetIndex());
672  CbmTrdPoint* point = (CbmTrdPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
673  id = point->GetTrackID();
674  //if (np > 1) cout << ip << " " << id << endl;
675  if (ids.find(id) == ids.end()) ids.insert(pair<Int_t,Int_t>(id,1));
676  else ++ids[id];
677  }
678  //if (np > 1) { cout << " digi " << digiM->GetNofLinks() << " " << ista << " " << hit->GetX() << " " << hit->GetY() << endl; exit(0); }
679  }
680  }
681  Int_t maxim = -1, idmax = -9;
682  for (map<Int_t,Int_t>::iterator it = ids.begin(); it != ids.end(); ++it) {
683  if (it->second > maxim) {
684  maxim = it->second;
685  idmax = it->first;
686  }
687  }
688  // Set vector ID as its flag
689  vec->SetFlag(idmax);
690  */
691 }
692 // -------------------------------------------------------------------------
693 
694 ClassImp(CbmTrdToTofVector);
CbmTrdToTofVector::fHitIds
std::map< Int_t, std::set< Int_t > > fHitIds
Definition: CbmTrdToTofVector.h:70
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmTrack::GetChiSq
Double_t GetChiSq() const
Definition: CbmTrack.h:58
CbmMatch
Definition: CbmMatch.h:22
CbmMCDataManager.h
CbmTrack::GetParamLast
const FairTrackParam * GetParamLast() const
Definition: CbmTrack.h:62
CbmTrdToTofVector::GetHits
void GetHits()
Definition: CbmTrdToTofVector.cxx:142
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmTrack.h:53
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmTrack::GetFlag
Int_t GetFlag() const
Definition: CbmTrack.h:57
CbmTrdToTofVector::fVectors
std::vector< CbmMuchTrack * > fVectors
Definition: CbmTrdToTofVector.h:71
CbmTrdToTofVector::fTrdTracks
TClonesArray * fTrdTracks
Definition: CbmTrdToTofVector.h:60
ClassImp
ClassImp(CbmTrdToTofVector)
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmTrdToTofVector::MatchTofToTrd
void MatchTofToTrd()
Definition: CbmTrdToTofVector.cxx:227
CbmTrdToTofVector::CbmTrdToTofVector
CbmTrdToTofVector()
Definition: CbmTrdToTofVector.cxx:36
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmMCDataArray.h
CbmTrack::SetPreviousTrackId
void SetPreviousTrackId(Int_t previousTrackId)
Definition: CbmTrack.h:72
CbmTrdToTofVector::RemoveClones
void RemoveClones()
Definition: CbmTrdToTofVector.cxx:377
CbmTrdToTofVector::FindLine
void FindLine(Int_t patt, Double_t *pars)
Definition: CbmTrdToTofVector.cxx:559
CbmTofDigi.h
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
CbmTrack::SetParamLast
void SetParamLast(const FairTrackParam *par)
Definition: CbmTrack.h:76
CbmTrack::SetNDF
void SetNDF(Int_t ndf)
Definition: CbmTrack.h:71
CbmMatch.h
CbmTrdToTofVector::SetTrackId
void SetTrackId(CbmMuchTrack *vec)
Definition: CbmTrdToTofVector.cxx:654
CbmMuchTrack
Definition: CbmMuchTrack.h:16
CbmTrdToTofVector::FindChi2
Double_t FindChi2(Int_t patt, Double_t *pars)
Definition: CbmTrdToTofVector.cxx:589
CbmTrdToTofVector::MakeVectors
void MakeVectors()
Definition: CbmTrdToTofVector.cxx:191
CbmTrack::SetFlag
void SetFlag(Int_t flag)
Definition: CbmTrack.h:69
CbmTrdToTofVector::AddVector
void AddVector(Int_t patt, Double_t chi2, Double_t *pars)
Definition: CbmTrdToTofVector.cxx:613
CbmTrdToTofVector::fHitX
std::multimap< Double_t, Int_t > fHitX
Definition: CbmTrdToTofVector.h:65
CbmTrdToTofVector::fZ
Double_t fZ[2]
Definition: CbmTrdToTofVector.h:61
CbmMuchTrack.h
CbmSetup.h
CbmTrack::SetMatch
void SetMatch(CbmMatch *match)
Definition: CbmTrack.cxx:80
CbmSetup::Instance
static CbmSetup * Instance()
Definition: CbmSetup.cxx:115
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmTrdToTofVector::fHitMatches
TClonesArray * fHitMatches
Definition: CbmTrdToTofVector.h:55
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmMatch::AddLink
void AddLink(const CbmLink &newLink)
Definition: CbmMatch.cxx:42
CbmTofAddress.h
CbmTrdToTofVector::Finish
virtual void Finish()
Definition: CbmTrdToTofVector.cxx:138
CbmTrack::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTrack.h:70
CbmSetup::GetGeoTag
Bool_t GetGeoTag(ECbmModuleId moduleId, TString &tag)
Definition: CbmSetup.cxx:101
CbmTrdToTofVector::fgkPlanes
static const Int_t fgkPlanes
Definition: CbmTrdToTofVector.h:49
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmTrdToTofVector::fHitY
std::multimap< Double_t, Int_t > fHitY
Definition: CbmTrdToTofVector.h:66
CbmTrdToTofVector::fHits
TClonesArray * fHits
Definition: CbmTrdToTofVector.h:54
CbmTofDigi
Data class for expanded digital TOF information.
Definition: CbmTofDigi.h:38
CbmTrdToTofVector::fHitTime
std::map< Int_t, Double_t > fHitTime
Definition: CbmTrdToTofVector.h:68
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
CbmTrdToTofVector::Init
virtual InitStatus Init()
Definition: CbmTrdToTofVector.cxx:57
CbmTrdToTofVector::fLineX
std::multimap< Double_t, TLine > fLineX
Definition: CbmTrdToTofVector.h:72
CbmTrdToTofVector::fPoints
CbmMCDataArray * fPoints
Definition: CbmTrdToTofVector.h:57
shape
UInt_t shape
Definition: CbmMvdSensorDigiToHitTask.cxx:73
CbmTrdToTofVector::~CbmTrdToTofVector
virtual ~CbmTrdToTofVector()
Definition: CbmTrdToTofVector.cxx:53
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofPoint.h
CbmTrdToTofVector::StoreVectors
void StoreVectors()
Definition: CbmTrdToTofVector.cxx:455
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTrdToTofVector::fTrackArray
TClonesArray * fTrackArray
Definition: CbmTrdToTofVector.h:52
CbmTrdToTofVector::fDigiMatches
TClonesArray * fDigiMatches
Definition: CbmTrdToTofVector.h:59
CbmTofPoint
Geometric intersection of a MC track with a TOFb detector.
Definition: CbmTofPoint.h:40
CbmTrdToTofVector::SetParContainers
virtual void SetParContainers()
Definition: CbmTrdToTofVector.cxx:100
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
CbmTrdToTofVector.h
CbmTrdToTofVector::Exec
virtual void Exec(Option_t *opt)
Definition: CbmTrdToTofVector.cxx:110
CbmTrdToTofVector
Definition: CbmTrdToTofVector.h:24
CbmTrdToTofVector::fDigis
TClonesArray * fDigis
Definition: CbmTrdToTofVector.h:58
CbmTrdToTofVector::CheckParams
void CheckParams()
Definition: CbmTrdToTofVector.cxx:344
CbmTofHit::GetPlaneId
Int_t GetPlaneId() const
Inherited from CbmBaseHit.
Definition: core/data/tof/CbmTofHit.h:88
CbmTrdToTofVector::ProcessPlane
void ProcessPlane(Int_t lay, Int_t patt, Int_t flag)
Definition: CbmTrdToTofVector.cxx:485
CbmDefs.h