CbmRoot
CbmTrdFindVectors.cxx
Go to the documentation of this file.
1 
5 #include "CbmTrdFindVectors.h"
6 #include "CbmDefs.h"
7 #include "CbmSetup.h"
8 //#include "CbmTrdModule.h"
9 //#include "CbmTrdDigiPar.h"
10 #include "CbmMCDataArray.h"
11 #include "CbmMCDataManager.h"
12 #include "CbmMatch.h"
13 #include "CbmMuchTrack.h"
14 #include "CbmTrdAddress.h"
15 #include "CbmTrdCluster.h"
16 #include "CbmTrdHit.h"
17 #include "CbmTrdPoint.h"
18 
19 //#include "FairRunAna.h"
20 //#include "FairRuntimeDb.h"
21 #include "FairRootManager.h"
22 
23 #include <TClonesArray.h>
24 #include <TGeoManager.h>
25 #include <TMath.h>
26 
27 #include <iostream>
28 
29 using std::cout;
30 using std::endl;
31 using std::map;
32 using std::multimap;
33 using std::pair;
34 
35 // ----- Default constructor -------------------------------------------
36 CbmTrdFindVectors::CbmTrdFindVectors()
37  : FairTask("TrdFindVectors")
38  , fTrackArray(NULL)
39  , fNofTracks(0)
40  , fHits(NULL)
41  , fClusters(NULL)
42  , fPoints(NULL)
43  , fDigiMatches(NULL)
44  , fZ0(99999)
45  , fErrX(1.0)
46  , fErrY(1.0)
47  , fCutChi2(24.0) // chi2/ndf = 6 for 4 hits
48 {}
49 // -------------------------------------------------------------------------
50 
51 // ----- Destructor ----------------------------------------------------
52 CbmTrdFindVectors::~CbmTrdFindVectors() {}
53 // -------------------------------------------------------------------------
54 
55 // ----- Public method Init (abstract in base class) --------------------
56 InitStatus CbmTrdFindVectors::Init() {
57 
58  // Get and check FairRootManager
59  FairRootManager* ioman = FairRootManager::Instance();
60  if (ioman == NULL) Fatal("Init", "RootManager not instantiated!");
61 
62  // Create and register TrdTrack array (if necessary)
63  //fTrackArray = static_cast<TClonesArray*> (ioman->GetObject("MuchVector"));
64  //if (fTrackArray == NULL) {
65  fTrackArray = new TClonesArray("CbmMuchTrack", 100);
66  ioman->Register("TrdVector", "Trd", fTrackArray, kTRUE);
67  //}
68  fHits = static_cast<TClonesArray*>(ioman->GetObject("TrdHit"));
69  fClusters = static_cast<TClonesArray*>(ioman->GetObject("TrdCluster"));
70  //fPoints = static_cast<TClonesArray*> (ioman->GetObject("TrdPoint"));
71  CbmMCDataManager* mcManager =
72  (CbmMCDataManager*) ioman->GetObject("MCDataManager");
73  if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
74  fPoints = mcManager->InitBranch("TrdPoint");
75  fDigiMatches = static_cast<TClonesArray*>(ioman->GetObject("TrdDigiMatch"));
76 
77  TString tag;
78  CbmSetup::Instance()->GetGeoTag(kTrd, tag);
79  cout << " ******* digiPar ******** " << tag << endl;
80  tag.Prepend("trd_");
81 
82  //CbmTrdGeoHandler geoHandler;
83  //geoHandler.Init();
84 
85  // Go to TRD
86  TGeoVolume* trdV = gGeoManager->GetVolume(tag);
87  TGeoNode* trd = NULL;
88  gGeoManager->CdTop();
89  TGeoNode* cave = gGeoManager->GetCurrentNode();
90  for (Int_t iNode = 0; iNode < cave->GetNdaughters(); iNode++) {
91  TGeoNode* node = cave->GetDaughter(iNode);
92  TString name = node->GetName();
93  if (name.Contains("trd")) {
94  trd = node;
95  gGeoManager->CdDown(iNode);
96  break;
97  }
98  }
99  if (!trd) {
100  cout << "-E- CbmTrdFindVectors::Init: Cannot find top TRD node" << endl;
101  return kFATAL;
102  }
103 
104  Int_t nlays = trdV->GetNdaughters();
105 
106  // Loop over layers
107  for (Int_t i = 0; i < nlays; ++i) {
108  TGeoVolume* lay = trdV->GetNode(i)->GetVolume();
109  if (!(TString(lay->GetName()).Contains("layer"))) continue;
110  gGeoManager->CdDown(i);
111  Int_t nmod = lay->GetNdaughters();
112 
113  // Loop over modules
114  for (Int_t j = 0; j < nmod; ++j) {
115  TGeoVolume* mod = lay->GetNode(j)->GetVolume();
116  if (!(TString(mod->GetName()).Contains("module"))) continue;
117  gGeoManager->CdDown(j);
118  Int_t nparts = mod->GetNdaughters();
119 
120  // Loop over module parts
121  for (Int_t k = 0; k < nparts; ++k) {
122  TGeoNode* part = mod->GetNode(k);
123  if (!(TString(part->GetName()).Contains("gas"))) continue;
124  gGeoManager->CdDown(k);
125  Double_t posLoc[3] = {0}, posGlob[3] = {0};
126  gGeoManager->LocalToMaster(posLoc, posGlob);
127  cout << " gas " << gGeoManager->GetPath() << " " << std::setprecision(4)
128  << posGlob[2] << endl;
129  if (i == 0) fZ0 = posGlob[2];
130  fDz[i] = posGlob[2] - fZ0;
131  gGeoManager->CdUp();
132  break;
133  }
134  gGeoManager->CdUp();
135  break;
136  } // for (Int_t j = 0; j < nmod;
137 
138  gGeoManager->CdUp();
139  } // for (Int_t i = 0; i < nlays;
140 
141  //exit(0);
142  ComputeMatrix(); // compute system matrices
143 
144  return kSUCCESS;
145 }
146 // -------------------------------------------------------------------------
147 
148 // ----- SetParContainers -------------------------------------------------
149 void CbmTrdFindVectors::SetParContainers() {
150  /*
151  fDigiPar = (CbmTrdDigiPar*) FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdDigiPar");
152  cout << " ******* digiPar ******** " << fDigiPar << endl;
153  exit(0);
154  */
155 }
156 // -------------------------------------------------------------------------
157 
158 // ----- Public method Exec --------------------------------------------
159 void CbmTrdFindVectors::Exec(Option_t* opt) {
160 
161  fTrackArray->Delete();
162 
163  Int_t nVecs = fVectors.size();
164  for (Int_t j = 0; j < nVecs; ++j)
165  delete fVectors[j];
166  fVectors.clear();
167 
168  // Do all processing
169 
170  // Get hits
171  GetHits();
172 
173  // Build vectors
174  MakeVectors();
175 
176  // Remove vectors with wrong orientation
177  // (using empirical cuts for omega muons at 8 GeV)
178  CheckParams();
179 
180  // Match vectors from 2 stations
181  //MatchVectors();
182 
183  // Remove clones
184  RemoveClones();
185 
186  // Store vectors
187  StoreVectors();
188 }
189 // -------------------------------------------------------------------------
190 
191 // ----- Public method Finish ------------------------------------------
192 void CbmTrdFindVectors::Finish() {}
193 // -------------------------------------------------------------------------
194 
195 // ----- Private method ComputeMatrix ----------------------------------
196 void CbmTrdFindVectors::ComputeMatrix() {
197  // Compute system matrices for different hit plane patterns
198 
199  Double_t dz2[fgkPlanes];
200  Bool_t onoff[fgkPlanes];
201 
202  for (Int_t i = 0; i < fgkPlanes; ++i) {
203  dz2[i] = fDz[i] * fDz[i];
204  onoff[i] = kTRUE;
205  }
206 
207  TMatrixD coef(4, 4);
208  Int_t pattMax = 1 << fgkPlanes, nDouble = 0, nTot = 0;
209 
210  // Loop over doublet patterns
211  for (Int_t ipat = 1; ipat < pattMax; ++ipat) {
212 
213  // Check if the pattern is valid: all doublets are active
214  nDouble = 0;
215  //for (Int_t j = 0; j < fgkPlanes; j += 2) if (ipat & (3 << j)) ++nDouble; else break;
216  for (Int_t j = 0; j < fgkPlanes; j += 1)
217  if (ipat & (1 << j)) ++nDouble;
218  if (nDouble < 3) continue;
219  ++nTot;
220 
221  for (Int_t j = 0; j < fgkPlanes; ++j)
222  onoff[j] = (ipat & (1 << j));
223 
224  coef = 0.0;
225  for (Int_t i = 0; i < fgkPlanes; ++i) {
226  if (!onoff[i]) continue;
227  coef(0, 0) += 1;
228  coef(0, 1) += 0;
229  coef(0, 2) += fDz[i];
230  coef(0, 3) += 0;
231  }
232  for (Int_t i = 0; i < fgkPlanes; ++i) {
233  if (!onoff[i]) continue;
234  coef(1, 0) += 0;
235  coef(1, 1) += 1;
236  coef(1, 2) += 0;
237  coef(1, 3) += fDz[i];
238  }
239  for (Int_t i = 0; i < fgkPlanes; ++i) {
240  if (!onoff[i]) continue;
241  coef(2, 0) += fDz[i];
242  coef(2, 1) += 0;
243  coef(2, 2) += dz2[i];
244  coef(2, 3) += 0;
245  }
246  for (Int_t i = 0; i < fgkPlanes; ++i) {
247  if (!onoff[i]) continue;
248  coef(3, 0) += 0;
249  coef(3, 1) += fDz[i];
250  coef(3, 2) += 0;
251  coef(3, 3) += dz2[i];
252  }
253 
254  TDecompLU* lu = new TDecompLU(4);
255  lu->SetMatrix(coef);
256  lu->SetTol(0.1 * lu->GetTol());
257  lu->Decompose();
258  //fLus.insert(pair<Int_t,TDecompLU*>(ipat,lu));
259  fLus[ipat] = lu;
260  TMatrixDSym cov(4);
261  cov.SetMatrixArray(coef.GetMatrixArray());
262  cov.Invert(); // covar. matrix
263  //fMatr.insert(pair<Int_t,TMatrixDSym*>(ipat,new TMatrixDSym(cov)));
264  fMatr[ipat] = new TMatrixDSym(cov);
265  TString buf = "";
266  for (Int_t jp = 0; jp < fgkPlanes; ++jp)
267  buf += Bool_t(ipat & (1 << jp));
268  Info("ComputeMatrix",
269  " Determinant: %s %i %f",
270  buf.Data(),
271  ipat,
272  coef.Determinant());
273  //if (ipat == 63) { coef.Print(); Info("ComputeMatrix"," Number of configurations: %i", nTot); }
274  //cout << " Number of configurations: " << nTot << endl; }
275  cov *= (1.2 * 1.2);
276  cout << TMath::Sqrt(cov(0, 0)) << " " << TMath::Sqrt(cov(1, 1)) << " "
277  << TMath::Sqrt(cov(2, 2)) << " " << TMath::Sqrt(cov(3, 3)) << endl;
278  cov *= (1.7 * 1.7 / 1.2 / 1.2);
279  cout << TMath::Sqrt(cov(0, 0)) << " " << TMath::Sqrt(cov(1, 1)) << " "
280  << TMath::Sqrt(cov(2, 2)) << " " << TMath::Sqrt(cov(3, 3)) << endl;
281  }
282 }
283 // -------------------------------------------------------------------------
284 
285 // ----- Private method GetHits ----------------------------------------
286 void CbmTrdFindVectors::GetHits() {
287  // Group hits according to their plane number
288 
289  for (Int_t i = 0; i < fgkPlanes; ++i) {
290  fHitPl[i].clear();
291  fHitX[i].clear();
292  }
293 
294  Int_t nHits = fHits->GetEntriesFast(), sel = 0;
295 
296  for (Int_t i = 0; i < nHits; ++i) {
297  CbmTrdHit* hit = (CbmTrdHit*) fHits->UncheckedAt(i);
298 
300  sel = 1; //SelectHitId(hit);
301  if (!sel) continue;
302  //
303  /*
304  Int_t address = hit->GetAddress();
305  Int_t layer = CbmTrdAddress::GetLayerId(address); // layer
306  Int_t colId = CbmTrdAddress::GetColumnId(address);
307  Int_t rowId = CbmTrdAddress::GetRowId(address);
308  */
309  Int_t layer = hit->GetPlaneId(), colId = 0;
310  fHitPl[layer].insert(pair<Int_t, Int_t>(colId, i));
311  fHitX[layer].insert(pair<Double_t, Int_t>(hit->GetX(), i));
312  }
313 }
314 // -------------------------------------------------------------------------
315 
316 // ----- Private method MakeVectors ------------------------------------
317 void CbmTrdFindVectors::MakeVectors() {
318  // Make vectors
319 
320  Int_t nvec = fVectors.size();
321  for (Int_t j = 0; j < nvec; ++j)
322  delete fVectors[j];
323  fVectors.clear();
324  CbmTrdHit* hit = NULL;
325 
326  cout << " TRD hits: " << fHits->GetEntriesFast() << endl;
327  for (Int_t lay = 0; lay < 2; ++lay) {
328  Int_t patt = 0, flag = 1, nhits = fHitPl[lay].size();
329  cout << " Hits: " << lay << " " << nhits << endl;
330  multimap<Int_t, Int_t>::iterator mit;
331 
332  for (mit = fHitPl[lay].begin(); mit != fHitPl[lay].end(); ++mit) {
333  Int_t indx = mit->second;
334  hit = (CbmTrdHit*) fHits->UncheckedAt(indx);
335  fXy[lay][0] = hit->GetX();
336  fXy[lay][1] = hit->GetY();
337  fXy[lay][2] = hit->GetDx();
338  fXy[lay][3] = hit->GetDy();
339  //fXy[lay][2] = TMath::Max (hit->GetDx(),fErrX[ista]);
340  //fXy[lay][3] = TMath::Max (hit->GetDy(),fErrY[ista]);
341  fXy[lay][4] = hit->GetZ();
342  //fXyi[lay][0] = CbmMuchAddress::GetModuleIndex(hit->GetAddress()); // sector No.
343  fXyi[lay][1] = indx;
344  patt = (1 << lay);
345  //cout << lay << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << " " << endl;
346  ProcessPlane(lay + 1, patt, flag);
347  }
348  }
349 }
350 // -------------------------------------------------------------------------
351 
352 // ----- Private method CheckParams ------------------------------------
353 void CbmTrdFindVectors::CheckParams() {
354  // Remove vectors with wrong orientation
355  // using empirical cuts for omega muons at 8 Gev
356 
357  const Double_t cut[2] = {0.6, 0.5}; // !!! empirical !!!
358 
359  Int_t nvec = fVectors.size();
360 
361  for (Int_t iv = 0; iv < nvec; ++iv) {
362  CbmMuchTrack* vec = fVectors[iv];
363  const FairTrackParam* params = vec->GetParamFirst();
364  Double_t dTx = params->GetTx() - params->GetX() / params->GetZ();
365  if (TMath::Abs(dTx) > cut[0])
366  vec->SetChiSq(-1.0);
367  else {
368  Double_t dTy = params->GetTy() - params->GetY() / params->GetZ();
369  if (TMath::Abs(dTy) > cut[1]) vec->SetChiSq(-1.0);
370  }
371  }
372 
373  for (Int_t iv = nvec - 1; iv >= 0; --iv) {
374  CbmMuchTrack* vec = fVectors[iv];
375  if (vec->GetChiSq() < 0) {
376  delete fVectors[iv];
377  fVectors.erase(fVectors.begin() + iv);
378  }
379  }
380  cout << " Vectors after parameter check: " << nvec << " " << fVectors.size()
381  << endl;
382 }
383 // -------------------------------------------------------------------------
384 
385 // ----- Private method RemoveClones -----------------------------------
386 void CbmTrdFindVectors::RemoveClones() {
387  // Remove clone vectors (having at least 3 the same hits)
388 
389  //Int_t nthr = 3, planes[20];
390  Int_t nthr = 2, planes[20];
391 
392  Int_t nvec = fVectors.size();
393 
394  // Do sorting according to "quality"
395  multimap<Double_t, CbmMuchTrack*> qMap;
396  multimap<Double_t, CbmMuchTrack*>::iterator it, it1;
397 
398  for (Int_t i = 0; i < nvec; ++i) {
399  CbmMuchTrack* vec = fVectors[i];
400  Double_t qual =
401  vec->GetNofHits() + (99 - TMath::Min(vec->GetChiSq(), 99.0)) / 100;
402  qMap.insert(pair<Double_t, CbmMuchTrack*>(-qual, vec));
403  }
404 
405  for (it = qMap.begin(); it != qMap.end(); ++it) {
406  CbmMuchTrack* vec = it->second;
407  if (vec->GetChiSq() < 0) continue;
408  for (Int_t j = 0; j < fgkPlanes; ++j)
409  planes[j] = -1;
410 
411  Int_t nhits = vec->GetNofHits();
412  for (Int_t ih = 0; ih < nhits; ++ih) {
413  CbmTrdHit* hit = (CbmTrdHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
414  Int_t lay = hit->GetPlaneId();
415  planes[lay] = vec->GetHitIndex(ih);
416  }
417 
418  it1 = it;
419  for (++it1; it1 != qMap.end(); ++it1) {
420  CbmMuchTrack* vec1 = it1->second;
421  if (vec1->GetChiSq() < 0) continue;
422  Int_t nsame = 0, same[fgkPlanes] = {0};
423 
424  Int_t nhits1 = vec1->GetNofHits();
425  //nthr = TMath::Min(nhits,nhits1) / 2;
426  //nthr = TMath::Min(nhits,nhits1) * 0.75;
427  for (Int_t ih = 0; ih < nhits1; ++ih) {
428  CbmTrdHit* hit = (CbmTrdHit*) fHits->UncheckedAt(vec1->GetHitIndex(ih));
429  Int_t lay = hit->GetPlaneId();
430  if (planes[lay] >= 0) {
431  if (vec1->GetHitIndex(ih) == planes[lay]) same[lay] = 1;
432  //else same[lay] = 0;
433  }
434  }
435  for (Int_t lay = 0; lay < fgkPlanes; ++lay)
436  nsame += same[lay];
437 
438  if (nsame >= nthr) {
439  // Too many the same hits
440  Int_t clone = 0;
441  if (nhits > nhits1 + 0)
442  clone = 1;
443  else if (vec->GetChiSq() * 1 <= vec1->GetChiSq())
444  clone = 1; // the same number of hits on 2 tracks
445  //else if (vec->GetChiSq() * 1.5 <= vec1->GetChiSq()) clone = 1; // the same number of hits on 2 tracks
446  if (clone) vec1->SetChiSq(-1.0);
447  }
448  }
449  } // for (it = qMap.begin();
450 
451  for (Int_t iv = nvec - 1; iv >= 0; --iv) {
452  CbmMuchTrack* vec = fVectors[iv];
453  if (vec->GetChiSq() < 0) {
454  delete fVectors[iv];
455  fVectors.erase(fVectors.begin() + iv);
456  }
457  }
458  cout << " Vectors after clones removed: " << nvec << " " << fVectors.size()
459  << endl;
460 }
461 // -------------------------------------------------------------------------
462 
463 // ----- Private method StoreVectors -----------------------------------
464 void CbmTrdFindVectors::StoreVectors() {
465  // Store vectors (CbmMuchTracks) into TClonesArray
466 
467  Int_t ntrs = fTrackArray->GetEntriesFast();
468  Int_t nHitsTot = fHits->GetEntriesFast();
469 
470  //set<Int_t> usedHits;
471  Int_t nvec = fVectors.size();
472 
473  for (Int_t iv = 0; iv < nvec; ++iv) {
474  CbmMuchTrack* tr =
475  new ((*fTrackArray)[ntrs++]) CbmMuchTrack(*(fVectors[iv]));
476  //cout << " Track: " << tr->GetNofHits() << endl;
477  //for (Int_t j = 0; j < tr->GetNofHits(); ++j) cout << j << " " << tr->GetHitIndex(j) << " " << fVectors[ist][iv]->GetHitIndex(j) << endl;
478  // Set hit flag (to check Lit tracking)
479  Int_t nhits = tr->GetNofHits();
480  }
481 }
482 // -------------------------------------------------------------------------
483 
484 // ----- Private method ProcessPlane -----------------------------------
485 void CbmTrdFindVectors::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.6,
489  0.5}; // !!! empirical !!! - cut on Tx(y) - x(y)/z
490  const Double_t cut1[2] = {
491  0.25, 0.15}; // !!! empirical !!! - cut on Tx(y) - Tx(y) (between segments)
492 
493  Double_t pars[4] = {0.0};
494  Int_t flag = 0;
495  multimap<Int_t, Int_t>::iterator mit;
496  // Clear bits
497  patt &= ~(1 << lay);
498 
499  for (mit = fHitPl[lay].begin(); mit != fHitPl[lay].end(); ++mit) {
500  Int_t indx = mit->second;
501 
502  CbmTrdHit* hit = (CbmTrdHit*) fHits->UncheckedAt(indx);
503  fXy[lay][0] = hit->GetX();
504  fXy[lay][1] = hit->GetY();
505  fXy[lay][2] = hit->GetDx();
506  fXy[lay][3] = hit->GetDy();
507  //fXy[lay][2] = TMath::Max (hit->GetDx(),fErrX[ista]);
508  //fXy[lay][3] = TMath::Max (hit->GetDy(),fErrY[ista]);
509  fXy[lay][4] = hit->GetZ();
510  // Clear bits
511  patt &= ~(1 << lay);
512 
513  // Check slopes
514  Int_t lay0 = 0;
515  if (!(patt & (1 << lay0))) ++lay0;
516  ;
517  if (TString::Itoa(patt, 2).CountChar('1') < 2) {
518  // First segment
519  Double_t dx = fXy[lay][0] - fXy[lay0][0];
520  Double_t dz = fXy[lay][4] - fXy[lay0][4];
521  //Double_t dTx = TMath::Abs(dx/dz) - TMath::Abs(fXy[lay][0]/fXy[lay][4]);
522  Double_t dTx = dx / dz - fXy[lay][0] / fXy[lay][4];
523  if (TMath::Abs(dTx) > cut[0]) continue;
524  Double_t dy = fXy[lay][1] - fXy[lay0][1];
525  //Double_t dTy = TMath::Abs(dy/dz) - TMath::Abs(fXy[lay][1]/fXy[lay][4]);
526  Double_t dTy = dy / dz - fXy[lay][1] / fXy[lay][4];
527  if (TMath::Abs(dTy) > cut[1]) continue;
528  } else {
529  Double_t dx = fXy[lay][0] - fXy[lay0][0];
530  Double_t dz = fXy[lay][4] - fXy[lay0][4];
531  Int_t lay1 = lay - 1;
532  if (!(patt & (1 << lay1))) --lay1;
533  ;
534  Double_t dx1 = fXy[lay][0] - fXy[lay1][0];
535  Double_t dz1 = fXy[lay][4] - fXy[lay1][4];
536  Double_t dTx = dx / dz - dx1 / dz1;
537  if (TMath::Abs(dTx) > cut1[0]) continue;
538  Double_t dy = fXy[lay][1] - fXy[lay0][1];
539  Double_t dy1 = fXy[lay][1] - fXy[lay1][1];
540  Double_t dTy = dy / dz - dy1 / dz1;
541  if (TMath::Abs(dTy) > cut1[1]) continue;
542  }
543 
544  fXyi[lay][1] = indx;
545  flag = 1;
546 
547  // Set bit
548  patt |= (1 << lay);
549  //cout << lay << " " << patt << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << " " << endl;
550 
551  if (lay < fgkPlanes - 1)
552  ProcessPlane(lay + 1, patt, flag);
553  else if (TString::Itoa(patt, 2).CountChar('1') > 2) {
554  // Straight line fit of the vector with > 2 hits
555  FindLine(patt, pars);
556  Double_t chi2 = FindChi2(patt, pars);
557  cout << " *** Stat: " << chi2 << " " << pars[0] << " " << pars[1] << endl;
558  if (chi2 <= fCutChi2)
559  AddVector(patt, chi2, pars); // add vector to the temporary container
560  }
561  } // for (mit = fHitPl[lay].begin();
562 
563  //if (flag == 0) {
564  if (flag >= 0) {
565  // Missing plane
566  if (flag0 == 0) return; // 2 consecutive missing planes
567  // Clear bits
568  patt &= ~(1 << lay);
569  if (lay < fgkPlanes - 1)
570  ProcessPlane(lay + 1, patt, flag);
571  else if (TString::Itoa(patt, 2).CountChar('1') > 2) {
572  // Straight line fit of the vector with > 2 hits
573  FindLine(patt, pars);
574  Double_t chi2 = FindChi2(patt, pars);
575  cout << " *** Stat: " << chi2 << " " << pars[0] << " " << pars[1] << endl;
576  if (chi2 <= fCutChi2)
577  AddVector(patt, chi2, pars); // add vector to the temporary container
578  }
579  }
580 }
581 // -------------------------------------------------------------------------
582 
583 // ----- Private method FindLine ---------------------------------------
584 void CbmTrdFindVectors::FindLine(Int_t patt, Double_t* pars) {
585  // Fit of hits to the straight line
586 
587  // Solve system of linear equations
588  Bool_t ok = kFALSE, onoff;
589  TVectorD b(4);
590  for (Int_t i = 0; i < fgkPlanes; ++i) {
591  onoff = patt & (1 << i);
592  if (!onoff) continue;
593  b[0] += fXy[i][0];
594  b[1] += fXy[i][1];
595  b[2] += fXy[i][0] * fDz[i];
596  b[3] += fXy[i][1] * fDz[i];
597  }
598 
599  //lu.Print();
600  TVectorD solve = fLus[patt]->Solve(b, ok);
601  //TVectorD solve = lu.TransSolve(b, ok);
602  //lu.Print();
603  for (Int_t i = 0; i < 4; ++i)
604  pars[i] = solve[i];
605  //for (Int_t i = 0; i < 4; ++i) { cout << pars[i] << " "; if (i == 3) cout << endl; }
606  //Double_t y1 = cosa / sina * (uz[1][0] * cosa - uz[0][0]) + uz[1][0] * sina;
607  //Double_t y2 = -cosa / sina * (uz[2][0] * cosa - uz[0][0]) - uz[2][0] * sina;
608  //cout << " OK " << ok << " " << y1 << " " << y2 << endl;
609 }
610 // -------------------------------------------------------------------------
611 
612 // ----- Private method FindChi2 ---------------------------------------
613 Double_t CbmTrdFindVectors::FindChi2(Int_t patt, Double_t* pars) {
614  // Compute chi2 of the fit
615 
616  Double_t chi2 = 0, x = 0, y = 0;
617  Bool_t onoff;
618 
619  for (Int_t i = 0; i < fgkPlanes; ++i) {
620  onoff = patt & (1 << i);
621  if (!onoff) continue;
622  x = pars[0] + pars[2] * fDz[i];
623  y = pars[1] + pars[3] * fDz[i];
624  //Double_t dx = (x - fXy[i][0]) / fErrX;
625  //Double_t dy = (y - fXy[i][1]) / fErrY;
626  Double_t dx = (x - fXy[i][0]) / TMath::Max(fXy[i][2], fErrX);
627  Double_t dy = (y - fXy[i][1]) / TMath::Max(fXy[i][3], fErrY);
628  chi2 += dx * dx;
629  chi2 += dy * dy;
630  }
631  //cout << " Chi2 = " << chi2 << endl;
632  return chi2;
633 }
634 // -------------------------------------------------------------------------
635 
636 // ----- Private method AddVector --------------------------------------
637 void CbmTrdFindVectors::AddVector(Int_t patt, Double_t chi2, Double_t* pars) {
638  // Add vector to the temporary container (as a MuchTrack)
639 
640  Bool_t refit = kFALSE; //kTRUE;
641  TMatrixDSym cov(4);
642 
643  if (refit) {
644  // Refit line with individual hit errors
645  //cout << " Before: " << chi2 << endl;
646  //Refit(patt, chi2, pars, cov);
647  //cout << " After: " << chi2 << endl;
648  } else {
649  cov = *fMatr[patt];
650  cov *= (fErrX * fErrX); // the same error in X and Y
651  }
652  cov.ResizeTo(5, 5);
653  cov(4, 4) = 1.0;
654 
655  CbmMuchTrack* track = new CbmMuchTrack();
656  track->SetChiSq(chi2);
657  FairTrackParam par(pars[0], pars[1], fZ0, pars[2], pars[3], 0.0, cov);
658  track->SetParamFirst(&par);
659  par.SetZ(fZ0 + fDz[fgkPlanes - 1]);
660  par.SetX(pars[0] + fDz[fgkPlanes - 1] * pars[2]);
661  par.SetY(pars[1] + fDz[fgkPlanes - 1] * pars[3]);
662  track->SetParamLast(&par);
663  //track->SetUniqueID(ista); // just to pass the value
664 
665  for (Int_t ipl = 0; ipl < fgkPlanes; ++ipl) {
666  if (!(patt & (1 << ipl))) continue;
667  track->AddHit(fXyi[ipl][1], kTRDHIT);
668  }
669  Int_t ndf = (track->GetNofHits() > 2) ? track->GetNofHits() * 2 - 4 : 1;
670  track->SetNDF(ndf);
671  SetTrackId(track); // set track ID as its flag
672  fVectors.push_back(track);
673 }
674 // -------------------------------------------------------------------------
675 
676 // ----- Private method SetTrackId -------------------------------------
677 void CbmTrdFindVectors::SetTrackId(CbmMuchTrack* vec) {
678  // Set vector ID as its flag (maximum track ID of its poins)
679 
680  map<Int_t, Int_t> ids;
681  Int_t nhits = vec->GetNofHits(), id = 0;
682 
683  //*
684  for (Int_t ih = 0; ih < nhits; ++ih) {
685  CbmTrdHit* hit = (CbmTrdHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
686  CbmTrdCluster* clus =
687  (CbmTrdCluster*) fClusters->UncheckedAt(hit->GetRefId());
688  Int_t nDigis = clus->GetNofDigis();
689 
690  for (Int_t j = 0; j < nDigis; ++j) {
691  CbmMatch* digiM = (CbmMatch*) fDigiMatches->UncheckedAt(clus->GetDigi(j));
692  Int_t np = digiM->GetNofLinks();
693 
694  for (Int_t ip = 0; ip < np; ++ip) {
695  CbmLink link = digiM->GetLink(ip);
696  //CbmTrdPoint* point = (CbmTrdPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
697  CbmTrdPoint* point = (CbmTrdPoint*) fPoints->Get(
698  link.GetFile(), link.GetEntry(), link.GetIndex());
699  id = point->GetTrackID();
700  //if (np > 1) cout << ip << " " << id << endl;
701  if (ids.find(id) == ids.end())
702  ids.insert(pair<Int_t, Int_t>(id, 1));
703  else
704  ++ids[id];
705  }
706  //if (np > 1) { cout << " digi " << digiM->GetNofLinks() << " " << ista << " " << hit->GetX() << " " << hit->GetY() << endl; exit(0); }
707  }
708  }
709  Int_t maxim = -1, idmax = -9;
710  for (map<Int_t, Int_t>::iterator it = ids.begin(); it != ids.end(); ++it) {
711  if (it->second > maxim) {
712  maxim = it->second;
713  idmax = it->first;
714  }
715  }
716  // Set vector ID as its flag
717  vec->SetFlag(idmax);
718  //*/
719 }
720 // -------------------------------------------------------------------------
721 
722 ClassImp(CbmTrdFindVectors);
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmTrdFindVectors::CheckParams
void CheckParams()
Definition: CbmTrdFindVectors.cxx:353
CbmTrack::GetChiSq
Double_t GetChiSq() const
Definition: CbmTrack.h:58
CbmMatch
Definition: CbmMatch.h:22
CbmMCDataManager.h
CbmTrdFindVectors::~CbmTrdFindVectors
virtual ~CbmTrdFindVectors()
Definition: CbmTrdFindVectors.cxx:52
CbmTrdFindVectors::fCutChi2
Double_t fCutChi2
Definition: CbmTrdFindVectors.h:127
CbmTrdFindVectors::FindLine
void FindLine(Int_t patt, Double_t *pars)
Definition: CbmTrdFindVectors.cxx:584
CbmTrdFindVectors::fDigiMatches
TClonesArray * fDigiMatches
Definition: CbmTrdFindVectors.h:116
CbmTrdAddress.h
Helper class to convert unique channel ID back and forth.
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
CbmTrdFindVectors::Exec
virtual void Exec(Option_t *opt)
Definition: CbmTrdFindVectors.cxx:159
CbmTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmTrack.h:53
ClassImp
ClassImp(CbmTrdFindVectors)
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmTrdFindVectors.h
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmTrdFindVectors::GetHits
void GetHits()
Definition: CbmTrdFindVectors.cxx:286
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmTrdHit
data class for a reconstructed Energy-4D measurement in the TRD
Definition: CbmTrdHit.h:35
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmTrdFindVectors::fMatr
std::map< Int_t, TMatrixDSym * > fMatr
Definition: CbmTrdFindVectors.h:129
CbmMCDataArray.h
CbmTrdFindVectors::RemoveClones
void RemoveClones()
Definition: CbmTrdFindVectors.cxx:386
CbmTrdFindVectors::fTrackArray
TClonesArray * fTrackArray
Definition: CbmTrdFindVectors.h:112
CbmTrdFindVectors::fZ0
Double_t fZ0
Definition: CbmTrdFindVectors.h:124
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
CbmTrack::SetParamLast
void SetParamLast(const FairTrackParam *par)
Definition: CbmTrack.h:76
CbmTrack::SetNDF
void SetNDF(Int_t ndf)
Definition: CbmTrack.h:71
CbmTrdFindVectors::SetParContainers
virtual void SetParContainers()
Definition: CbmTrdFindVectors.cxx:149
CbmMatch.h
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmMuchTrack
Definition: CbmMuchTrack.h:16
CbmTrack::SetFlag
void SetFlag(Int_t flag)
Definition: CbmTrack.h:69
CbmCluster::GetNofDigis
Int_t GetNofDigis() const
Number of digis in cluster.
Definition: CbmCluster.h:69
CbmTrdFindVectors::fXyi
Double_t fXyi[fgkPlanes][3]
Definition: CbmTrdFindVectors.h:120
CbmTrdCluster
Data Container for TRD clusters.
Definition: CbmTrdCluster.h:23
CbmTrdFindVectors::CbmTrdFindVectors
CbmTrdFindVectors()
Definition: CbmTrdFindVectors.cxx:36
CbmMuchTrack.h
CbmSetup.h
CbmTrdFindVectors::FindChi2
Double_t FindChi2(Int_t patt, Double_t *pars)
Definition: CbmTrdFindVectors.cxx:613
CbmTrdFindVectors::Finish
virtual void Finish()
Definition: CbmTrdFindVectors.cxx:192
CbmTrdFindVectors::Init
virtual InitStatus Init()
Definition: CbmTrdFindVectors.cxx:56
CbmSetup::Instance
static CbmSetup * Instance()
Definition: CbmSetup.cxx:115
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmTrdFindVectors::ComputeMatrix
void ComputeMatrix()
Definition: CbmTrdFindVectors.cxx:196
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmTrdFindVectors::fErrX
Double_t fErrX
Definition: CbmTrdFindVectors.h:125
CbmTrdFindVectors::fErrY
Double_t fErrY
Definition: CbmTrdFindVectors.h:126
CbmTrack::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTrack.h:70
CbmTrdHit.h
Class for hits in TRD detector.
CbmTrdFindVectors::fHits
TClonesArray * fHits
Definition: CbmTrdFindVectors.h:114
CbmSetup::GetGeoTag
Bool_t GetGeoTag(ECbmModuleId moduleId, TString &tag)
Definition: CbmSetup.cxx:101
CbmTrdFindVectors::ProcessPlane
void ProcessPlane(Int_t lay, Int_t patt, Int_t flag)
Definition: CbmTrdFindVectors.cxx:485
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
kTRDHIT
@ kTRDHIT
Definition: CbmHit.h:25
CbmTrdPoint
Definition: CbmTrdPoint.h:23
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
ECbmModuleId::kTrd
@ kTrd
Transition Radiation Detector.
CbmTrdFindVectors::fClusters
TClonesArray * fClusters
Definition: CbmTrdFindVectors.h:115
CbmTrdFindVectors::fPoints
CbmMCDataArray * fPoints
Definition: CbmTrdFindVectors.h:118
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTrdFindVectors::MakeVectors
void MakeVectors()
Definition: CbmTrdFindVectors.cxx:317
CbmTrdFindVectors::AddVector
void AddVector(Int_t patt, Double_t chi2, Double_t *pars)
Definition: CbmTrdFindVectors.cxx:637
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTrdPoint.h
CbmTrack::SetParamFirst
void SetParamFirst(const FairTrackParam *par)
Definition: CbmTrack.h:75
CbmTrdFindVectors::fDz
Double_t fDz[fgkPlanes]
Definition: CbmTrdFindVectors.h:121
CbmTrdFindVectors::fVectors
std::vector< CbmMuchTrack * > fVectors
Definition: CbmTrdFindVectors.h:134
CbmTrdCluster.h
Data Container for TRD clusters.
CbmTrdFindVectors::fgkPlanes
static const Int_t fgkPlanes
Definition: CbmTrdFindVectors.h:56
CbmTrdFindVectors::StoreVectors
void StoreVectors()
Definition: CbmTrdFindVectors.cxx:464
CbmTrdFindVectors::fLus
std::map< Int_t, TDecompLU * > fLus
Definition: CbmTrdFindVectors.h:123
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
CbmTrdFindVectors::fHitX
std::multimap< Double_t, Int_t > fHitX[fgkPlanes]
Definition: CbmTrdFindVectors.h:133
CbmTrdFindVectors::SetTrackId
void SetTrackId(CbmMuchTrack *vec)
Definition: CbmTrdFindVectors.cxx:677
CbmTrdFindVectors
Definition: CbmTrdFindVectors.h:29
CbmTrack::AddHit
void AddHit(Int_t index, HitType type)
Definition: CbmTrack.cxx:75
CbmTrdFindVectors::fXy
Double_t fXy[fgkPlanes][5]
Definition: CbmTrdFindVectors.h:119
CbmTrdFindVectors::fHitPl
std::multimap< Int_t, Int_t > fHitPl[fgkPlanes]
Definition: CbmTrdFindVectors.h:131
CbmDefs.h
CbmTrdHit::GetPlaneId
Int_t GetPlaneId() const
Inherited from CbmBaseHit.
Definition: CbmTrdHit.h:73