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