CbmRoot
CbmMuchMergeVectors.cxx
Go to the documentation of this file.
1 
5 #include "CbmMuchMergeVectors.h"
6 #include "CbmKFTrack.h"
7 #include "CbmMuchDigiMatch.h"
8 #include "CbmMuchFindVectorsGem.h"
9 #include "CbmMuchModule.h"
10 #include "CbmMuchPoint.h"
11 #include "CbmMuchStation.h"
12 #include "CbmMuchTrack.h"
13 #include "CbmStsTrack.h"
14 //#include "CbmTrackMatch.h"
15 #include "CbmTrackMatchNew.h"
16 
17 #include "FairEventHeader.h"
18 #include "FairField.h"
19 #include "FairLogger.h"
20 #include "FairRootManager.h"
21 
22 #include <TClonesArray.h>
23 #include <TGeoManager.h>
24 #include <TMath.h>
25 #include <TMatrixD.h>
26 #include <TMatrixFLazy.h>
27 #include <TROOT.h>
28 #include <TVectorD.h>
29 
30 #include <iostream>
31 
32 using std::cout;
33 using std::endl;
34 using std::map;
35 using std::multimap;
36 using std::pair;
37 using std::set;
38 
39 FILE* lun1 = 0x0; //fopen("chi21.dat","w");
40 
41 // ----- Default constructor -------------------------------------------
42 CbmMuchMergeVectors::CbmMuchMergeVectors()
43  : FairTask("MuchMergeVectors")
44  , fGeoScheme(CbmMuchGeoScheme::Instance())
45  , fTrackArray(NULL)
46  , fNofTracks(0)
47  , fHits(NULL)
48  , fGemHits(NULL)
49  , fPoints(NULL)
50  , fDigiMatches(NULL)
51  , fVecArray(NULL)
52  , fTracksSts(NULL)
53  , fTrStsMatch(NULL)
54  , fTracksLit(NULL)
55  , fStatFirst(-1)
56 //fCutChi2(40)//200.0)
57 {
58  for (Int_t i = 0; i < 9; ++i)
59  fCutChi2[i] = 40;
60  fCutChi2[0] *= 2;
61  //fCutChi2[fgkStat-2] *= 2;
62 }
63 // -------------------------------------------------------------------------
64 
65 // ----- Destructor ----------------------------------------------------
66 CbmMuchMergeVectors::~CbmMuchMergeVectors() {
67  fTrackArray->Delete();
68 
69  for (map<Int_t, TMatrixDSym*>::iterator it = fMatr.begin(); it != fMatr.end();
70  ++it)
71  delete it->second;
72 }
73 // -------------------------------------------------------------------------
74 
75 // ----- Public method Init (abstract in base class) --------------------
76 InitStatus CbmMuchMergeVectors::Init() {
77 
78  // Get and check FairRootManager
79  FairRootManager* ioman = FairRootManager::Instance();
80  if (ioman == NULL) Fatal("Init", "RootManager not instantiated!");
81 
82  // Create and register MuchTrack array
83  fTrackArray = new TClonesArray("CbmMuchTrack", 100);
84  ioman->Register("MuchVectorTrack", "Much", fTrackArray, kTRUE);
85 
86  CbmMuchFindHitsStraws* hitFinder =
87  (CbmMuchFindHitsStraws*) FairRun::Instance()->GetTask(
88  "CbmMuchFindHitsStraws");
89  //if (hitFinder == NULL) Fatal("Init", "CbmMuchFindHitsStraws not run!");
90  //fDiam = hitFinder->GetDiam(0);
91 
92  fGemHits = static_cast<TClonesArray*>(ioman->GetObject("MuchPixelHit"));
93  fPoints = static_cast<TClonesArray*>(ioman->GetObject("MuchPoint"));
94  fVecArray = static_cast<TClonesArray*>(ioman->GetObject("MuchVector"));
95  fTracksSts = static_cast<TClonesArray*>(ioman->GetObject("StsTrack"));
96  fTrStsMatch = static_cast<TClonesArray*>(ioman->GetObject("StsTrackMatch"));
97  fTracksLit = static_cast<TClonesArray*>(ioman->GetObject("MuchTrack"));
98 
99  // Find first straw station and get some geo constants
100  Int_t nSt = fGeoScheme->GetNStations();
101  for (Int_t i = 0; i < nSt; ++i) {
102  CbmMuchStation* st = fGeoScheme->GetStation(i);
103  CbmMuchModule* mod = fGeoScheme->GetModule(i, 0, 0, 0);
104  fRmin[i] = st->GetRmin();
105  fRmax[i] = st->GetRmax();
106  if (mod->GetDetectorType() == 2) {
107  if (fStatFirst < 0)
108  fStatFirst = CbmMuchAddress::GetStationIndex(st->GetDetectorId());
109  //cout << " First station: " << fStatFirst << endl;
110  Int_t nLays = st->GetNLayers();
111  fNdoubl = nLays;
112  for (Int_t lay = 0; lay < nLays; ++lay) {
113  CbmMuchLayer* layer = st->GetLayer(lay);
114  Double_t phi = hitFinder->GetPhi(lay);
115  for (Int_t iside = 0; iside < 2; ++iside) {
116  CbmMuchLayerSide* side = layer->GetSide(iside);
117  Int_t plane = lay * 2 + iside;
118  fDz[plane] = side->GetZ();
119  fCosa[plane] = TMath::Cos(phi);
120  fSina[plane] = TMath::Sin(phi);
121  if (plane == 0) fZ0[i - fStatFirst] = side->GetZ();
122  }
123  }
124  }
125  }
126  for (Int_t i = fgkPlanes - 1; i >= 0; --i) {
127  fDz[i] -= fDz[0];
128  //cout << fDz[i] << " ";
129  }
130  //cout << endl;
131 
132  // Get absorbers from GEM vector finder
133  FairTask* vecFinder = FairRun::Instance()->GetTask("VectorFinder");
134  vecFinder->GetListOfTasks()->ls();
135  CbmMuchFindVectorsGem* gemFinder =
136  (CbmMuchFindVectorsGem*) vecFinder->GetListOfTasks()->FindObject(
137  "MuchFindVectorsGem");
138  if (gemFinder == NULL) Fatal("Init", "CbmMuchFindVectorsGem not run!");
139 
140  Int_t nAbs = gemFinder->GetAbsorbers(fZabs0, fX0abs);
141 
142  cout << " \n !!! MUCH Absorbers: " << nAbs << "\n Zbeg, Zend, X0:";
143  for (Int_t j = 0; j < nAbs; ++j)
144  cout << " " << std::setprecision(4) << fZabs0[j][0] << ", " << fZabs0[j][1]
145  << ", " << fX0abs[j] << ";";
146  cout << endl << endl;
147 
148  return kSUCCESS;
149 }
150 // -------------------------------------------------------------------------
151 
152 // ----- SetParContainers -------------------------------------------------
153 void CbmMuchMergeVectors::SetParContainers() {}
154 // -------------------------------------------------------------------------
155 
156 // ----- Public method Exec --------------------------------------------
157 void CbmMuchMergeVectors::Exec(Option_t* opt) {
158 
159  //gErrorIgnoreLevel = kInfo; //debug level
160  //gErrorIgnoreLevel = kWarning; //debug level
161  gLogger->SetLogScreenLevel("INFO");
162  //gLogger->SetLogScreenLevel("WARNING");
163 
164  fTrackArray->Delete();
165 
166  // Do all processing
167 
168  // Get vectors
169  GetVectors();
170 
171  // Match vectors
172  MatchVectors();
173 
174  // Merge vectors
175  MergeVectors();
176 
177  // Select final tracks - remove ghosts
178  SelectTracks();
179 
180  // Add first station (as a filter)
181  //AddStation1();
182 }
183 // -------------------------------------------------------------------------
184 
185 // ----- Public method Finish ------------------------------------------
186 void CbmMuchMergeVectors::Finish() {
187  fTrackArray->Clear();
188 
189  for (map<Int_t, TMatrixDSym*>::iterator it = fMatr.begin(); it != fMatr.end();
190  ++it)
191  delete it->second;
192 }
193 // -------------------------------------------------------------------------
194 
195 // ----- Private method GetVectors -------------------------------------
196 void CbmMuchMergeVectors::GetVectors() {
197  // Group vectors according to the station number
198 
199  std::map<Int_t, CbmMuchTrack*>::iterator it;
200  for (Int_t ista = 0; ista < fgkStat; ++ista) {
201  if (ista == 0) {
202  // STS tracks
203  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it)
204  if (it->first >= 0) delete it->second;
205  }
206  fVectors[ista].clear();
207  }
208 
209  Int_t nVecs = fVecArray->GetEntriesFast(), sel = 0;
210  for (Int_t i = 0; i < nVecs; ++i) {
211  CbmMuchTrack* vec = (CbmMuchTrack*) fVecArray->UncheckedAt(i);
212 
214  sel = 1; //SelectHitId(hit);
215  if (!sel) continue;
216  //
217 
218  Int_t ista = vec->GetUniqueID() + 1; // offset due to STS tracks
219  fVectors[ista].insert(pair<Int_t, CbmMuchTrack*>(i, vec));
220  }
221 
222  // Get STS tracks and extrapolate them to the first absorber face
223  if (fTracksSts == NULL) return;
224  Int_t nSts = fTracksSts->GetEntriesFast(), ista = 0;
225  FairTrackParam param, parRear;
226  Double_t pMin = 0.25;
227  //Double_t pMin = 0.7; //AZ
228 
229  for (Int_t i = 0; i < nSts; ++i) {
230  CbmStsTrack* tr = (CbmStsTrack*) fTracksSts->UncheckedAt(i);
231  if (tr->GetNofHits() < 4) continue; // too few hits
232  Double_t ppp = 1. / TMath::Abs(tr->GetParamLast()->GetQp());
233  //if (ppp < pMin + 0.1) continue; // too low momentum
234  if (ppp < pMin + 0.05) continue; // too low momentum
235  CbmKFTrack kfTrack = CbmKFTrack(*tr, 0);
236  if (kfTrack.Extrapolate(fZabs0[ista][0]))
237  continue; // extrapolation error to absorber front face
238  kfTrack.GetTrackParam(param);
239  if (kfTrack.Extrapolate(fZabs0[ista][1]))
240  continue; // extrapolation error to absorber rear face
241  kfTrack.GetTrackParam(parRear);
242  Double_t r2 = param.GetX() * param.GetX() + param.GetY() * param.GetY();
243  if (TMath::Sqrt(r2) > fRmax[0]) continue; // outside outer acceptance
244  // Take into account track angle
245  Double_t cos2th =
246  1.0 + param.GetTx() * param.GetTx() + param.GetTy() * param.GetTy();
247  //if (ppp / TMath::Sqrt(cos2th) < pMin + 0.1) continue; // too low momentum
248  if (ppp / TMath::Sqrt(cos2th) < pMin + 0.05) continue; // too low momentum
249  CbmMuchTrack* vec = new CbmMuchTrack();
250  param.SetQp(1.0 / ppp);
251  vec->SetParamFirst(&param);
252  vec->SetPreviousTrackId(i); // store STS track index
253  vec->SetUniqueID(9); // station "No. 9"
254  parRear.SetQp(TMath::Abs(parRear.GetQp()));
255  vec->SetParamLast(&parRear);
256  //CbmTrackMatch *trMatch = (CbmTrackMatch*) fTrStsMatch->UncheckedAt(i);
257  //vec->SetFlag(trMatch->GetMCTrackId());
258  if (fTrStsMatch) {
259  CbmTrackMatchNew* trMatch =
260  (CbmTrackMatchNew*) fTrStsMatch->UncheckedAt(i);
261  vec->SetFlag(trMatch->GetMatchedLink().GetIndex());
262  //} else vec->SetFlag(SetStsTrackId(i));
263  } else if (tr->GetMatch() && tr->GetMatch()->GetNofLinks()) {
264  vec->SetFlag(tr->GetMatch()->GetMatchedLink().GetIndex());
265  } else
266  vec->SetFlag(-1);
267  fVectors[ista].insert(pair<Int_t, CbmMuchTrack*>(i + nVecs, vec));
268  //cout << " Before absorber: " << i << " " << nVecs << " " << i+nVecs << endl;
269  }
270 
271  // Extrapolate STS tracks through the first absorber
272  TMatrixF matr = TMatrixF(5, 5);
273  TMatrixF unit(TMatrixF::kUnit, matr);
274  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it) {
275  CbmMuchTrack* tr = it->second;
276  FairTrackParam parFirst = *tr->GetParamFirst();
277  FairTrackParam parLast = *tr->GetParamLast();
278  parFirst = parLast;
279  TMatrixFSym cov(5);
280  parFirst.CovMatrix(cov);
281  /*
282  Double_t zbeg = parFirst.GetZ();
283  Double_t dz = fZabs0[0][1] - zbeg;
284  // Propagate params
285  //parFirst.SetX(parFirst.GetX() + dz * parFirst.GetTx());
286  //parFirst.SetY(parFirst.GetY() + dz * parFirst.GetTy());
287  //parFirst.SetZ(parFirst.GetZ() + dz);
288  parFirst.SetX(parLast.GetX());
289  parFirst.SetY(parLast.GetY());
290  parFirst.SetZ(parLast.GetZ());
291  parFirst.SetTx(parLast.GetTx());
292  parFirst.SetTy(parLast.GetTy());
293  TMatrixFSym cov(5);
294  parFirst.CovMatrix(cov);
295  cov(4,4) = 1.0;
296  TMatrixF ff = unit;
297  ff(2,0) = ff(3,1) = dz;
298  TMatrixF cf(cov,TMatrixF::kMult,ff);
299  TMatrixF fcf(ff,TMatrixF::kTransposeMult,cf);
300  cov.SetMatrixArray(fcf.GetMatrixArray());
301  */
302  PassAbsorber(0, fZabs0[0], fX0abs[0], parFirst, cov, 1);
303  parFirst.SetCovMatrix(cov);
304  //tr->SetParamLast(&parFirst);
305  tr->SetParamFirst(&parFirst);
306  //cout << " After absorber: " << it->first << " " << parFirst.GetX() << " " << parFirst.GetY() << endl;
307  }
308 }
309 // -------------------------------------------------------------------------
310 
311 // ----- Private method SetStsTrackId ----------------------------------
312 /*
313 Int_t CbmMuchMergeVectors::SetStsTrackId(Int_t indx)
314 {
315  // Set STS track ID
316 
317  CbmStsTrack *tr = (CbmStsTrack*) fTracksSts->UncheckedAt(indx);
318  Int_t nHits = tr->GetNofStsHits();
319 
320  for (Int_t i = 0; i < nHits; ++i) {
321  CbmStsHit *hit = fStsHits->UncheckedAt(tr->GetStsHitIndex(i));
322  }
323 }
324 */
325 // -------------------------------------------------------------------------
326 
327 // ----- Private method MatchVectors -----------------------------------
328 void CbmMuchMergeVectors::MatchVectors() {
329  // Match vectors (CbmMuchTracks) from 2 stations going upstream -
330  // remove vectors without matching
331 
332  const Double_t window0 = 7.0; //10.0; //
333  TMatrixF matr = TMatrixF(5, 5);
334  TMatrixF unit(TMatrixF::kUnit, matr);
335  map<Int_t, CbmMuchTrack*>::iterator it, it1;
336  multimap<Double_t, pair<Int_t, CbmMuchTrack*>> xMap[2];
337  multimap<Double_t, pair<Int_t, CbmMuchTrack*>>::iterator mit, mit1, mitb,
338  mite;
339  map<CbmMuchTrack*, Int_t> matchOK[2];
340  TMatrixFSym cov(5);
341 
342  //AZ for (Int_t iabs = 2; iabs >= 0; --iabs) { // last absorber has been checked already
343  for (Int_t iabs = 3; iabs >= 0;
344  --iabs) { // last absorber has been checked already
345  Double_t window = window0;
346  if (iabs == fgkStat - 2) window = window0 * 2.0;
347  //Int_t ibeg = fTrackArray->GetEntriesFast();
348  xMap[0].clear();
349  xMap[1].clear();
350  matchOK[0].clear();
351  matchOK[1].clear();
352 
353  for (Int_t ist = 0; ist < 2; ++ist) {
354  // Propagate vectors to the absorber faces
355  Int_t ista = iabs - 1 + ist + 1;
356 
357  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it) {
358  CbmMuchTrack* tr = it->second;
359  FairTrackParam parFirst = *tr->GetParamFirst();
360  parFirst.CovMatrix(cov);
361  cov(4, 4) = 1.0;
362  if (ista > 0) { // STS tracks has been propagated already
363  Double_t zbeg = parFirst.GetZ();
364  Double_t dz = fZabs0[iabs][1] - zbeg;
365  // Propagate params
366  parFirst.SetX(parFirst.GetX() + dz * parFirst.GetTx());
367  parFirst.SetY(parFirst.GetY() + dz * parFirst.GetTy());
368  parFirst.SetZ(parFirst.GetZ() + dz);
369  TMatrixF ff = unit;
370  ff(2, 0) = ff(3, 1) = dz;
371  //cout << " Cov: " << dz << " " << cov(0,0) << " " << cov(1,1) << endl;
372  TMatrixF cf(cov, TMatrixF::kMult, ff);
373  TMatrixF fcf(ff, TMatrixF::kTransposeMult, cf);
374  cov.SetMatrixArray(fcf.GetMatrixArray());
375  if (ist == 0 && iabs > 0)
376  PassAbsorber(
377  ist + iabs * 2, fZabs0[iabs], fX0abs[iabs], parFirst, cov, 0);
378  //if (ist) cov.Print(); ;
379  }
380  cov.Invert(); // weight matrix
381  parFirst.SetCovMatrix(cov);
382  //tr->SetParamFirst(&parFirst);
383  tr->SetParamLast(&parFirst);
384  xMap[ist].insert(
385  pair<Double_t, pair<Int_t, CbmMuchTrack*>>(parFirst.GetX(), *it));
386  matchOK[ist][tr] = -1;
387  }
388  } // for (Int_t ist = 0; ist < 2;
389 
390  Int_t ista0 = iabs - 1 + 1, ista1 = iabs + 1;
391 
392  for (mit = xMap[0].begin(); mit != xMap[0].end(); ++mit) {
393  CbmMuchTrack* tr1 = mit->second.second;
394  FairTrackParam parOk = *tr1->GetParamLast();
395  FairTrackParam par1 = *tr1->GetParamLast();
396  TMatrixFSym w1(5);
397  par1.CovMatrix(w1);
398  Float_t pars1[5] = {(Float_t) par1.GetX(),
399  (Float_t) par1.GetY(),
400  (Float_t) par1.GetTx(),
401  (Float_t) par1.GetTy(),
402  1.0};
403  TMatrixF p1(5, 1, pars1);
404  TMatrixF wp1(w1, TMatrixF::kTransposeMult, p1);
405  Double_t x0 = parOk.GetX(), y0 = parOk.GetY();
406  mitb = xMap[1].lower_bound(x0 - window); // lower X-window edge
407  mite = xMap[1].upper_bound(x0 + window); // upper X-window edge
408 
409  for (mit1 = mitb; mit1 != mite; ++mit1) {
410  CbmMuchTrack* tr2 = mit1->second.second;
411  FairTrackParam par2 = *tr2->GetParamLast();
412  if (par2.GetY() - y0 < -window) continue;
413  if (par2.GetY() - y0 > window) continue;
414  TMatrixFSym w2(5);
415  par2.CovMatrix(w2);
416  TMatrixFSym w20 = w2;
417  Float_t pars2[5] = {(Float_t) par2.GetX(),
418  (Float_t) par2.GetY(),
419  (Float_t) par2.GetTx(),
420  (Float_t) par2.GetTy(),
421  1.0};
422  TMatrixF p2(5, 1, pars2);
423  TMatrixF wp2(w2, TMatrixF::kTransposeMult, p2);
424  wp2 += wp1;
425  w2 += w1;
426  w2.Invert();
427  TMatrixF p12(w2, TMatrixF::kTransposeMult, wp2);
428 
429  // Compute Chi2
430  TMatrixF p122 = p12;
431  TMatrixF pMerge = p12;
432  p12 -= p1;
433  TMatrixF wDp1(w1, TMatrixF::kMult, p12);
434  TMatrixF chi21(p12, TMatrixF::kTransposeMult, wDp1);
435  p122 -= p2;
436  TMatrixF wDp2(w20, TMatrixF::kMult, p122);
437  TMatrixF chi22(p122, TMatrixF::kTransposeMult, wDp2);
438  Double_t c2 = chi21(0, 0) + chi22(0, 0);
439  //cout << " Chi2: " << chi21(0,0) << " " << chi22(0,0) << " " << c2 << endl;
440  if (c2 < 0 || c2 > fCutChi2[iabs]) continue;
441  matchOK[0][tr1] = 1;
442  matchOK[1][tr2] = 1;
443  // Merged track parameters
444  /*
445  parOk.SetX(pMerge(0,0));
446  parOk.SetY(pMerge(1,0));
447  parOk.SetZ(par2.GetZ());
448  parOk.SetTx(pMerge(2,0));
449  parOk.SetTy(pMerge(3,0));
450  parOk.SetCovMatrix(w2);
451  */
452  //AddTrack(ista0, tr1, tr2, mit->second.first, mit1->second.first, parOk, c2); // add track
453  } // for (mit1 = xMap[1].begin();
454  } // for (mit = xMap[0].begin();
455 
456  for (Int_t ist = 0; ist < 2; ++ist) {
457  Int_t ista = iabs - 1 + ist + 1;
458 
459  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it) {
460  CbmMuchTrack* tr = it->second;
461  if (matchOK[ist][tr] > 0) continue;
462  if (ista == 0) delete tr;
463  fVectors[ista].erase(it);
464  }
465  cout << " MergeVectors: vectors after matching in station " << ista
466  << ": " << fVectors[ista].size() << endl;
467  }
468  } // for (Int_t iabs = 2; iabs >= 0;
469 }
470 // -------------------------------------------------------------------------
471 
472 // ----- Private method MergeVectors -----------------------------------
473 void CbmMuchMergeVectors::MergeVectors() {
474  // Match vectors (CbmMuchTracks) from 2 stations
475 
476  const Double_t window0 = 5.0; //10.0; //
477  TMatrixF matr = TMatrixF(5, 5);
478  TMatrixF unit(TMatrixF::kUnit, matr);
479  map<Int_t, CbmMuchTrack*>::iterator it, it1;
480  multimap<Double_t, pair<Int_t, CbmMuchTrack*>> xMap[2];
481  multimap<Double_t, pair<Int_t, CbmMuchTrack*>>::iterator mit, mit1;
482 
483  for (Int_t iabs = 0; iabs <= 3; ++iabs) {
484  Int_t ibeg = fTrackArray->GetEntriesFast();
485  xMap[0].clear();
486  xMap[1].clear();
487 
488  for (Int_t ist = 0; ist < 2; ++ist) {
489  // Propagate vectors to the absorber faces
490  Int_t ista = iabs - 1 + ist + 1;
491  //if (iabs == 1 && ist == 0) --ista; // !!! take extrapolated STS tracks - skip first station
492  Int_t imerged = (fVectors[ista].begin()->first < 0)
493  ? 1
494  : 0; // merged vectors stored with negative station No.
495  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it) {
496  if (imerged && it->first >= 0)
497  break; // for merged vectors: skip original (unmerged) vectors
498  CbmMuchTrack* tr = it->second;
499  FairTrackParam parFirst = *tr->GetParamFirst();
500  //if (ist == 0) parFirst = *tr->GetParamLast();
501  Double_t zbeg = parFirst.GetZ();
502  //Double_t dz = zAbs0[iabs][ist] - zbeg;
503  Double_t dz = fZabs0[iabs][1] - zbeg;
504  //Double_t dz = fZabs0[iabs][0] - zbeg;
505  //Double_t dz = (zAbs0[iabs][0] + zAbs0[iabs][1]) / 2 - zbeg;
506  // Propagate params
507  parFirst.SetX(parFirst.GetX() + dz * parFirst.GetTx());
508  parFirst.SetY(parFirst.GetY() + dz * parFirst.GetTy());
509  parFirst.SetZ(parFirst.GetZ() + dz);
510  TMatrixFSym cov(5);
511  parFirst.CovMatrix(cov);
512  cov(4, 4) = 1.0;
513  //cov.Print();
514  TMatrixF ff = unit;
515  //ff.Print();
516  //ff(0,2) = ff(1,3) = dz;
517  ff(2, 0) = ff(3, 1) = dz;
518  //cout << " Cov: " << dz << " " << cov(0,0) << " " << cov(1,1) << endl;
519  TMatrixF cf(cov, TMatrixF::kMult, ff);
520  TMatrixF fcf(ff, TMatrixF::kTransposeMult, cf);
521  cov.SetMatrixArray(fcf.GetMatrixArray());
522  //cout << " Cov1: " << dz << " " << cov(0,0) << " " << cov(1,1) << endl;
523  // Go thru absorber
524  //if (ist) { cov.Print(); cout << " -0- " << endl; }
525  //AZ if (ist == 0) {
526  if (ist == 0 && iabs > 0) {
527  if (iabs == 1)
528  PassAbsorber(ist + iabs * 2,
529  fZabs0[iabs],
530  fX0abs[iabs],
531  parFirst,
532  cov,
533  1); // STS track
534  else
535  PassAbsorber(
536  ist + iabs * 2, fZabs0[iabs], fX0abs[iabs], parFirst, cov, 1);
537  }
538  //if (ist) cov.Print(); ;
539  cov.Invert(); // weight matrix
540  parFirst.SetCovMatrix(cov);
541  //tr->SetParamFirst(&parFirst);
542  tr->SetParamLast(&parFirst);
543  xMap[ist].insert(
544  pair<Double_t, pair<Int_t, CbmMuchTrack*>>(parFirst.GetX(), *it));
545 
546  //Info("MergeVectors", "Absorber %i, station %i, trID %i, X = %f, Y = %f, Z = %f", iabs,ista,
547  // tr->GetFlag(),parFirst.GetX(),parFirst.GetY(),parFirst.GetZ());
548  //cov.Print();
549  }
550  } // for (Int_t ist = 0; ist < 2;
551 
552  Int_t ista0 = iabs - 1 + 1, ista1 = iabs + 1;
553  //if (iabs == 1 && ista0 == 1) --ista0; // !!! take extrapolated STS tracks - skip first station
554  map<Int_t, CbmMuchTrack*> mapMerged;
555  Double_t window = window0;
556  if (iabs == 3) window *= 2;
557 
558  for (mit = xMap[0].begin(); mit != xMap[0].end(); ++mit) {
559  //Int_t imerged0 = (fVectors[ista0].begin()->first < 0) ? 1: 0;
560  //if (iabs && !imerged0) break; //!!! no track merging occured for previous absorber
561  //Int_t imerged1 = (fVectors[ista1].begin()->first < 0) ? 1: 0;
562  //if (imerged0 && it->first >= 0) break; // for merged tracks: exclude original vectors
563  CbmMuchTrack* tr1 = mit->second.second;
564  FairTrackParam parOk = *tr1->GetParamLast();
565  //if (1./parOk.GetQp() < 0.1) continue; // too low momentum
566  if (1. / parOk.GetQp() < 0.05) continue; // too low momentum
567  FairTrackParam par1 = *tr1->GetParamLast();
568  TMatrixFSym w1(5);
569  par1.CovMatrix(w1);
570  Float_t pars1[5] = {(Float_t) par1.GetX(),
571  (Float_t) par1.GetY(),
572  (Float_t) par1.GetTx(),
573  (Float_t) par1.GetTy(),
574  1.0};
575  TMatrixF p1(5, 1, pars1);
576  TMatrixF wp1(w1, TMatrixF::kTransposeMult, p1);
577  Double_t x0 = parOk.GetX(), y0 = parOk.GetY();
578 
579  for (mit1 = xMap[1].begin(); mit1 != xMap[1].end(); ++mit1) {
580  //if (imerged1 && it1->first >= 0) break; // for merged tracks: exclude original vectors
581  CbmMuchTrack* tr2 = mit1->second.second;
582  //if (tr2->GetUniqueID() == 2 && fNdoubl == 3 && 1./parOk.GetQp() < 0.05+0.37) break; // exclude low-mom. tracks - for merged straw stations (0.37 - energy loss in last absorber)
583  //FairTrackParam par2 = *tr2->GetParamFirst();
584  FairTrackParam par2 = *tr2->GetParamLast();
585  if (par2.GetX() - x0 < -window) continue;
586  if (par2.GetX() - x0 > window) break;
587  if (par2.GetY() - y0 < -window) continue;
588  if (par2.GetY() - y0 > window) continue;
589  TMatrixFSym w2(5);
590  par2.CovMatrix(w2);
591  TMatrixFSym w20 = w2;
592  Float_t pars2[5] = {(Float_t) par2.GetX(),
593  (Float_t) par2.GetY(),
594  (Float_t) par2.GetTx(),
595  (Float_t) par2.GetTy(),
596  1.0};
597  TMatrixF p2(5, 1, pars2);
598  TMatrixF wp2(w2, TMatrixF::kTransposeMult, p2);
599  wp2 += wp1;
600  w2 += w1;
601  w2.Invert();
602  TMatrixF p12(w2, TMatrixF::kTransposeMult, wp2);
603  //p12.Print();
604 
605  // Compute Chi2
606  TMatrixF p122 = p12;
607  TMatrixF pMerge = p12;
608  p12 -= p1;
609  TMatrixF wDp1(w1, TMatrixF::kMult, p12);
610  TMatrixF chi21(p12, TMatrixF::kTransposeMult, wDp1);
611  p122 -= p2;
612  TMatrixF wDp2(w20, TMatrixF::kMult, p122);
613  TMatrixF chi22(p122, TMatrixF::kTransposeMult, wDp2);
614  Double_t c2 = chi21(0, 0) + chi22(0, 0);
615  //cout << " Chi2: " << chi21(0,0) << " " << chi22(0,0) << " " << c2 << endl;
616  if (c2 < 0 || c2 > fCutChi2[iabs]) continue;
617  // Merged track parameters
618  parOk.SetX(pMerge(0, 0));
619  parOk.SetY(pMerge(1, 0));
620  parOk.SetZ(par2.GetZ());
621  parOk.SetTx(pMerge(2, 0));
622  parOk.SetTy(pMerge(3, 0));
623  parOk.SetCovMatrix(w2);
624  AddTrack(ista0,
625  tr1,
626  tr2,
627  mit->second.first,
628  mit1->second.first,
629  parOk,
630  c2); // add track
631  Int_t evNo = FairRun::Instance()->GetEventHeader()->GetMCEntryNumber();
632  //fprintf(lun1,"%6d %6d %6d %10.3e \n",evNo,tr1->GetFlag(),tr2->GetFlag(),c2);
633  /* Debug
634  TMatrixFSym covA = w1;
635  covA.Invert();
636  if (lun1 && tr1->GetFlag() == tr2->GetFlag() && tr1->GetFlag() < 2) fprintf(lun1,"%6d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",evNo,par1.GetX()-par2.GetX(),par1.GetY()-par2.GetY(),par1.GetTx()-par2.GetTx(),par1.GetTy()-par2.GetTy(),TMath::Sqrt(covA(0,0)),TMath::Sqrt(covA(1,1)),TMath::Sqrt(covA(2,2)),TMath::Sqrt(covA(3,3)));
637  */
638  } // for (it1 = fVectors[ista1].begin();
639  } // for (it = fVectors[ista0].begin();
640 
641  // Remove clones for each absorber
642  //RemoveClones(ibeg, iabs, mapMerged);
643  RemoveClones(ibeg, -1, mapMerged);
644  //if (iabs == 3) RemoveClones(ibeg, -1, mapMerged); // !!! remove clones
645  //else RemoveClones(ibeg, 1, mapMerged); // !!! do not remove clones
646  if (mapMerged.size() == 0) break; // no merging with STS tracks !!!
647 
648  // Add merged tracks
649  for (it = mapMerged.begin(); it != mapMerged.end(); ++it) {
650  fVectors[ista1][it->first] = it->second;
651  }
652 
653  } // for (Int_t iabs = 0; iabs <= 3;
654 }
655 // -------------------------------------------------------------------------
656 
657 // ----- Private method PassAbsorber -----------------------------------
658 void CbmMuchMergeVectors::PassAbsorber(Int_t ist,
659  Double_t* zabs,
660  Double_t x0,
661  FairTrackParam& parFirst,
662  TMatrixFSym& cov,
663  Int_t pFlag) {
664  // Go thru absorber
665 
666  Double_t x = parFirst.GetX();
667  Double_t y = parFirst.GetY();
668  Double_t z = parFirst.GetZ();
669  Double_t tx = parFirst.GetTx();
670  Double_t ty = parFirst.GetTy();
671  //Double_t dz = zabs[0] - z;
672  Double_t dz = (zabs[0] - zabs[1]) / 1;
673 
674  // Propagate params
675  //parFirst.SetX(x + dz * tx);
676  //parFirst.SetY(y + dz * ty);
677  //parFirst.SetZ(z + dz);
678  //return;
679 
680  Double_t aaa = 1.0 + tx * tx + ty * ty;
681  dz *= TMath::Sqrt(aaa); //
682  Double_t l = TMath::Abs(dz);
683  Double_t l2 = l * l;
684  Double_t l3 = l2 * l;
685 
686  TMatrixFSym covMS(5);
687  //covMS(0,0) = l3 / 3 * (1.0 + tx * tx);
688  //covMS(1,1) = l3 / 3 * (1.0 + ty * ty);
689  //covMS(2,2) = l * (1.0 + tx * tx);
690  //covMS(3,3) = l * (1.0 + ty * ty);
691  covMS(0, 0) = l2 / 3 * (1.0 + tx * tx);
692  covMS(1, 1) = l2 / 3 * (1.0 + ty * ty);
693  covMS(2, 2) = 1 * (1.0 + tx * tx);
694  covMS(3, 3) = 1 * (1.0 + ty * ty);
695 
696  //covMS(0,1) = covMS(1,0) = l3 / 3 * tx * ty;
697  covMS(0, 1) = covMS(1, 0) = l2 / 3 * tx * ty;
698 
699  Int_t icor = (ist % 2 == 0) ? 1 : -1;
700  //icor = -1; //
701 
702  //covMS(0,2) = covMS(2,0) = icor * l2 / 2 * (1.0 + tx * tx);
703  //covMS(1,2) = covMS(2,1) = icor * l2 / 2 * tx * ty;
704  covMS(0, 2) = covMS(2, 0) = icor * l / 2 * (1.0 + tx * tx);
705  covMS(1, 2) = covMS(2, 1) = icor * l / 2 * tx * ty;
706 
707  //covMS(0,3) = covMS(3,0) = icor * l2 / 2 * tx * ty;
708  //covMS(1,3) = covMS(3,1) = icor * l2 / 2 * (1.0 + ty * ty);
709  //covMS(2,3) = covMS(3,2) = l * tx * ty;
710  covMS(0, 3) = covMS(3, 0) = icor * l / 2 * tx * ty;
711  covMS(1, 3) = covMS(3, 1) = icor * l / 2 * (1.0 + ty * ty);
712  covMS(2, 3) = covMS(3, 2) = 1 * tx * ty;
713 
714  covMS *= aaa;
715 
716  Double_t dxx0 = l / x0, angle = 0.0;
717  if (pFlag == 0)
718  angle =
719  0.0136 / 0.985 / 0.6
720  * (1. + 0.038 * TMath::Log(dxx0)); // mu at p = 0.6 - peak at 8A GeV
721  //if (pFlag == 0) angle = 0.0136 / 0.997 / 1.3 * (1. + 0.038 * TMath::Log(dxx0)); // mu at p = 1.3 GeV
722  else {
723  Double_t pmom = 1.0 / TMath::Abs(parFirst.GetQp());
724  Double_t beta = pmom / TMath::Sqrt(pmom * pmom + 0.106 * 0.106);
725  angle = 0.0136 / beta / pmom * (1. + 0.038 * TMath::Log(dxx0));
726  }
727  //covMS *= (angle * angle / x0);
728  covMS *= (angle * angle * dxx0);
729  cov += covMS;
730  if (pFlag == 0) return; // no momentum update
731 
732  // Momentum update due to energy loss
733  Double_t pLoss[6] = {0.25, 0.25, 0.25, 0.37, 0, 0};
734  Double_t ppp = 1. / parFirst.GetQp();
735  ppp -= (pLoss[ist / 2] * TMath::Sqrt(aaa)); // ppp is always non-negative
736  parFirst.SetQp(1. / ppp);
737 }
738 // -------------------------------------------------------------------------
739 
740 // ----- Private method AddTrack ---------------------------------------
741 void CbmMuchMergeVectors::AddTrack(Int_t ista0,
742  CbmMuchTrack* tr1,
743  CbmMuchTrack* tr2,
744  Int_t indx1,
745  Int_t indx2,
746  FairTrackParam& parOk,
747  Double_t c2) {
748  // Store merged vector (CbmMuchTracks) into TClonesArray
749 
750  Int_t ntrs = fTrackArray->GetEntriesFast();
751 
752  CbmMuchTrack* track = new ((*fTrackArray)[ntrs]) CbmMuchTrack();
753  //track->SetParamFirst(tr2->GetParamFirst());
754  track->SetParamFirst(&parOk);
755  if (tr1->GetUniqueID() == -9) {
756  // Exclude STS track chi2
757  track->SetChiSq(c2 + tr2->GetChiSq());
758  track->SetNDF(4 + tr2->GetNDF());
759  } else {
760  track->SetChiSq(c2 + tr1->GetChiSq() + tr2->GetChiSq());
761  track->SetNDF(4 + tr1->GetNDF() + tr2->GetNDF());
762  }
763  track->SetUniqueID(tr2->GetUniqueID()); // station number
764  //track->SetPreviousTrackId(tr1->GetPreviousTrackId());
765  track->SetPreviousTrackId(indx1);
766  if (indx1 >= 0) {
767  if (ista0) track->AddHit(indx1, kMUCHSTRAWHIT); // add vector index
768  //else track->AddHit(indx1, kHIT); // add STS vector index
769  else
770  track->AddHit(tr1->GetPreviousTrackId(), kHIT); // add STS track index
771  //if (indx2 < 0) track->SetPreviousTrackId(-indx2 - 1); // index of previous track
772  } else {
773  // Merged track
774  Int_t nmerged = tr1->GetNofHits();
775  for (Int_t j = 0; j < nmerged; ++j)
776  track->AddHit(tr1->GetHitIndex(j),
777  tr1->GetHitType(j)); // add vector index
778  if (indx1 < 0)
779  track->SetPreviousTrackId(-indx1 - 1); // index of previous track
780  }
781  if (indx2 >= 0)
782  track->AddHit(indx2, kMUCHSTRAWHIT); // add vector index
783  else {
784  // Merged track
785  Int_t nmerged = tr2->GetNofHits();
786  for (Int_t j = 0; j < nmerged; ++j)
787  track->AddHit(tr2->GetHitIndex(j),
788  tr2->GetHitType(j)); // add vector index
789  }
790  if (tr1->GetFlag() == tr2->GetFlag())
791  track->SetFlag(tr1->GetFlag());
792  else
793  track->SetFlag(-1);
794 
795  //Info("AddTrack", "trID1=%i, trID2=%i, chi2=%f, merged vectors %i", tr1->GetFlag(),tr2->GetFlag(),track->GetChiSq(),track->GetNofHits());
796  gLogger->Info(MESSAGE_ORIGIN,
797  "CbmMuchMergeVectors::AddTrack: ista=%i, trID1=%i, trID2=%i, "
798  "chi2=%f, merged vectors %i",
799  ista0,
800  tr1->GetFlag(),
801  tr2->GetFlag(),
802  track->GetChiSq(),
803  track->GetNofHits());
804 }
805 // -------------------------------------------------------------------------
806 
807 // ----- Private method RemoveClones -----------------------------------
808 void CbmMuchMergeVectors::RemoveClones(Int_t ibeg,
809  Int_t iabs,
810  map<Int_t, CbmMuchTrack*>& mapMerged) {
811  // Remove clone tracks (having at least 1 the same vector)
812 
813  Int_t ntrs = fTrackArray->GetEntriesFast();
814 
815  if (iabs != 1) { // Do not remove clones when merging stations 0 and 1
816  multimap<Double_t, Int_t> c2merge;
817  for (Int_t i1 = ibeg; i1 < ntrs; ++i1) {
818  CbmMuchTrack* tr1 = (CbmMuchTrack*) fTrackArray->UncheckedAt(i1);
819  /*
820  Double_t qual = 0.0;
821  if (tr1->GetUniqueID() == 9) {
822  CbmMuchTrack *trLit = (CbmMuchTrack*) fTracksLit->UncheckedAt(tr1->GetHitIndex(0));
823  qual = trLit->GetNofHits() + (499 - TMath::Min(tr1->GetChiSq(),499.0)) / 500;
824  } else qual = tr1->GetNofHits() + (499 - TMath::Min(tr1->GetChiSq(),499.0)) / 500;
825  */
826  Double_t qual =
827  tr1->GetNofHits()
828  + (499 - TMath::Min(tr1->GetChiSq() / tr1->GetNDF(), 499.0)) / 500;
829  c2merge.insert(pair<Double_t, Int_t>(-qual, i1));
830  }
831 
832  multimap<Double_t, Int_t>::iterator it, it1;
833  for (it = c2merge.begin(); it != c2merge.end(); ++it) {
834  CbmMuchTrack* tr1 = (CbmMuchTrack*) fTrackArray->UncheckedAt(it->second);
835  if (tr1 == NULL) continue;
836  Int_t nvecs1 = tr1->GetNofHits();
837 
838  it1 = it;
839  for (++it1; it1 != c2merge.end(); ++it1) {
840  CbmMuchTrack* tr2 =
841  (CbmMuchTrack*) fTrackArray->UncheckedAt(it1->second);
842  if (tr2 == NULL) continue;
843  Int_t nvecs2 = tr2->GetNofHits();
844  if (tr2->GetUniqueID() != tr1->GetUniqueID()) continue;
845 
846  Bool_t over = kFALSE;
847  for (Int_t iv1 = 0; iv1 < nvecs1; ++iv1) {
848  for (Int_t iv2 = iv1; iv2 < nvecs2; ++iv2) {
849  //for (Int_t iv1 = 0; iv1 < 2; ++iv1) {
850  //for (Int_t iv2 = iv1; iv2 < 2; ++iv2) {
851  if (iv2 != iv1) continue;
852  if (tr1->GetHitType(iv1) != tr2->GetHitType(iv2)) continue;
853  //if (tr1->GetUniqueID() == 9 && iv1 == 1) continue; // !!! share vectors in 1st station
854  if (tr1->GetHitIndex(iv1) != tr2->GetHitIndex(iv2)) continue;
855  // Count number of overlaps for the given vector
856  if (iv2) {
857  CbmMuchTrack* vec =
858  (CbmMuchTrack*) fVecArray->UncheckedAt(tr2->GetHitIndex(iv2));
859  Int_t clones = vec->GetPreviousTrackId();
860  ++clones;
861  if (clones == 0) ++clones;
862  vec->SetPreviousTrackId(clones);
863  }
864  gLogger->Info(MESSAGE_ORIGIN,
865  "CbmMuchMergeVectors:RemoveClones: qual1 %f, qual2 "
866  "%f, trID1 %i, trID2 %i, ista %i, p1 %f, p2 %f",
867  it->first,
868  it1->first,
869  tr1->GetFlag(),
870  tr2->GetFlag(),
871  iv1,
872  1 / tr1->GetParamFirst()->GetQp(),
873  1 / tr2->GetParamFirst()->GetQp());
874  fTrackArray->RemoveAt(it1->second);
875  over = kTRUE;
876  break;
877  }
878  if (over) break;
879  }
880  }
881  }
882  fTrackArray->Compress();
883  } // if (iabs != 1)
884 
885  ntrs = fTrackArray->GetEntriesFast();
886  // Add track to the map (with negative index)
887  for (Int_t i1 = ibeg; i1 < ntrs; ++i1) {
888  CbmMuchTrack* tr1 = (CbmMuchTrack*) fTrackArray->UncheckedAt(i1);
889  mapMerged[-i1 - 1] = tr1;
890  }
891 }
892 // -------------------------------------------------------------------------
893 
894 // ----- Private method SelectTracks -----------------------------------
895 void CbmMuchMergeVectors::SelectTracks() {
896  // Remove ghost tracks (having at least N the same hits (i.e. fired tubes))
897 
898  const Int_t nMax[2] = {2, 2}, nPl = 40;
899  //Int_t nVecMin = (fNdoubl == 3) ? 4 : 5;
900  Int_t nVecMin = 4;
901  Int_t planes[nPl], ntrs = fTrackArray->GetEntriesFast();
902 
903  multimap<Double_t, Int_t> c2merge;
904  for (Int_t i = 0; i < ntrs; ++i) {
905  CbmMuchTrack* tr = (CbmMuchTrack*) fTrackArray->UncheckedAt(i);
906  if (tr->GetNofHits() < nVecMin) continue;
907  /*
908  Double_t qual = 0.0;
909  if (tr->GetUniqueID() == 9) {
910  CbmMuchTrack *trLit = (CbmMuchTrack*) fTracksLit->UncheckedAt(tr->GetHitIndex(0));
911  qual = trLit->GetNofHits() + (499 - TMath::Min(tr->GetChiSq(),499.0)) / 500;
912  } else qual = tr->GetNofHits() + (499 - TMath::Min(tr->GetChiSq(),499.0)) / 500;
913  */
914  Double_t qual =
915  tr->GetNofHits()
916  + (499 - TMath::Min(tr->GetChiSq() / tr->GetNDF(), 499.0)) / 500;
917  c2merge.insert(pair<Double_t, Int_t>(-qual, i));
918  }
919  if (c2merge.size() < 2) return;
920 
921  multimap<Double_t, Int_t>::iterator it, it1;
922  for (it = c2merge.begin(); it != c2merge.end(); ++it) {
923  CbmMuchTrack* tr1 = (CbmMuchTrack*) fTrackArray->UncheckedAt(it->second);
924  if (tr1 == NULL) continue;
925  Int_t nvecs1 = tr1->GetNofHits();
926  for (Int_t j = 0; j < nPl; ++j)
927  planes[j] = -1;
928 
929  for (Int_t iv = 0; iv < nvecs1; ++iv) {
930  if (tr1->GetHitType(iv) == kHIT) {
931  // STS track
932  planes[nPl - 1] = tr1->GetHitIndex(iv);
933  continue;
934  }
935  CbmMuchTrack* vec =
936  (CbmMuchTrack*) fVecArray->UncheckedAt(tr1->GetHitIndex(iv));
937  TClonesArray* hits = fHits;
938  if (vec->GetUniqueID() < 2 || fStatFirst < 0) hits = fGemHits;
939  Int_t nh = vec->GetNofHits();
940  for (Int_t ih = 0; ih < nh; ++ih) {
941  CbmHit* hit = (CbmHit*) hits->UncheckedAt(vec->GetHitIndex(ih));
942  Int_t ipl = hit->GetPlaneId() - 1;
943  planes[ipl] = vec->GetHitIndex(ih);
944  }
945  }
946 
947  it1 = it;
948  for (++it1; it1 != c2merge.end(); ++it1) {
949  CbmMuchTrack* tr2 = (CbmMuchTrack*) fTrackArray->UncheckedAt(it1->second);
950  if (tr2 == NULL) continue;
951  //if (tr2->GetUniqueID() != tr1->GetUniqueID()) continue;
952  Int_t nvecs2 = tr2->GetNofHits(), nover[2] = {0};
953  Bool_t over = kFALSE;
954 
955  for (Int_t iv = 0; iv < nvecs2; ++iv) {
956  if (tr2->GetHitType(iv) == kHIT) {
957  // STS track
958  if (planes[nPl - 1] >= 0 && planes[nPl - 1] == tr2->GetHitIndex(iv)) {
959  // The same STS track
960  gLogger->Info(MESSAGE_ORIGIN,
961  "Track quality: qual1 %f, qual2 %f, trID1 %i, trID2 "
962  "%i, the same STS track: %i",
963  it->first,
964  it1->first,
965  tr1->GetFlag(),
966  tr2->GetFlag(),
967  tr2->GetHitIndex(iv));
968  fTrackArray->RemoveAt(it1->second);
969  break;
970  }
971  if (tr2->GetUniqueID() != tr1->GetUniqueID()) break;
972  continue;
973  }
974  CbmMuchTrack* vec =
975  (CbmMuchTrack*) fVecArray->UncheckedAt(tr2->GetHitIndex(iv));
976  Int_t nh = vec->GetNofHits();
977  TClonesArray* hits = fHits;
978  if (vec->GetUniqueID() < 2 || fStatFirst < 0) hits = fGemHits;
979  for (Int_t ih = 0; ih < nh; ++ih) {
980  CbmHit* hit = (CbmHit*) hits->UncheckedAt(vec->GetHitIndex(ih));
981  Int_t ipl = hit->GetPlaneId() - 1;
982  if (planes[ipl] < 0) continue;
983  if (planes[ipl] == vec->GetHitIndex(ih)) {
984  if (hits == fGemHits)
985  ++nover[0];
986  else
987  ++nover[1];
988  if (nover[0] >= nMax[0] || nover[1] >= nMax[1]) {
989  //cout << ipl << " " << vec->GetHitIndex(ih) << endl;
990  gLogger->Info(MESSAGE_ORIGIN,
991  "Track quality: qual1 %f, qual2 %f, trID1 %i, "
992  "trID2 %i, overlaps: %i, %i",
993  it->first,
994  it1->first,
995  tr1->GetFlag(),
996  tr2->GetFlag(),
997  nover[0],
998  nover[1]);
999  fTrackArray->RemoveAt(it1->second);
1000  over = kTRUE;
1001  break;
1002  }
1003  }
1004  }
1005  if (over) break;
1006  } // for (Int_t iv = 0; iv < nvecs2;
1007 
1008  } // for (++it1; it1 != c2merge.end();
1009  }
1010  fTrackArray->Compress();
1011 }
1012 // -------------------------------------------------------------------------
1013 
1014 // ------ Private method AddStation1 -----------------------------------
1015 void CbmMuchMergeVectors::AddStation1() {
1016  // Add vector from the first station as a filter
1017 
1018  const Int_t nVecsMin = 4;
1019  map<Int_t, CbmMuchTrack*>::iterator it, it1;
1020  TMatrixF matr = TMatrixF(5, 5);
1021  TMatrixF unit(TMatrixF::kUnit, matr);
1022 
1023  if (fVectors[0].size() == 0) return;
1024  //if (fVectors[0].begin()->first >= 0) return; // no merged tracks were found
1025  if (fVectors[0].begin()->second->GetNofHits() != nVecsMin)
1026  return; // no long tracks were found
1027 
1028  for (Int_t ist = 0; ist < 1; ++ist) {
1029 
1030  // Propagate vectors to the absorber faces
1031  Int_t iabs = 0, ista = 1;
1032  for (it = fVectors[ista].begin(); it != fVectors[ista].end(); ++it) {
1033  //if (imerged && it->first >= 0) break; // for merged vectors: skip original (unmerged) vectors
1034  CbmMuchTrack* tr = it->second;
1035  FairTrackParam parFirst = *tr->GetParamFirst();
1036  //if (ist == 0) parFirst = *tr->GetParamLast();
1037  Double_t zbeg = parFirst.GetZ();
1038  Double_t dz = fZabs0[iabs][1] - zbeg;
1039  // Propagate params
1040  parFirst.SetX(parFirst.GetX() + dz * parFirst.GetTx());
1041  parFirst.SetY(parFirst.GetY() + dz * parFirst.GetTy());
1042  parFirst.SetZ(parFirst.GetZ() + dz);
1043  TMatrixFSym cov(5);
1044  parFirst.CovMatrix(cov);
1045  cov(4, 4) = 1.0;
1046  TMatrixF ff = unit;
1047  ff(2, 0) = ff(3, 1) = dz;
1048  TMatrixF cf(cov, TMatrixF::kMult, ff);
1049  TMatrixF fcf(ff, TMatrixF::kTransposeMult, cf);
1050  cov.SetMatrixArray(fcf.GetMatrixArray());
1051  cov.Invert(); // weight matrix
1052  parFirst.SetCovMatrix(cov);
1053  tr->SetParamFirst(&parFirst);
1054  } // for (it = fVectors[ista].begin();
1055  } // for (Int_t ist = 0; ist < 1;
1056 
1057  // Get STS tracks extrapolated through the first absorber
1058  Int_t nvecs = fVecArray->GetEntriesFast();
1059  for (it = fVectors[0].begin(); it != fVectors[0].end(); ++it) {
1060  if (it->first >= 0) break; // no merged tracks anymore
1061  Int_t indSts = (it->second)->GetHitIndex(0); // STS track index
1062  CbmMuchTrack* tr1 = fVectors[0][nvecs + indSts];
1063  //cout << " STS index: " << indSts << " " << tr1->GetParamFirst()->GetX() << " " << tr1->GetParamFirst()->GetY() << endl;
1064 
1065  // Merge with vectors from first station
1066  FairTrackParam parOk = *tr1->GetParamFirst();
1067  FairTrackParam par1 = *tr1->GetParamFirst();
1068  TMatrixFSym w1(5);
1069  par1.CovMatrix(w1);
1070  w1.Invert(); // weight matrix
1071  Float_t pars1[5] = {(Float_t) par1.GetX(),
1072  (Float_t) par1.GetY(),
1073  (Float_t) par1.GetTx(),
1074  (Float_t) par1.GetTy(),
1075  1.0};
1076  TMatrixF p1(5, 1, pars1);
1077  TMatrixF wp1(w1, TMatrixF::kTransposeMult, p1);
1078  Double_t c2min = 999999.0;
1079  //cout << " STS " << endl;
1080  //par1.Print();
1081  //cout << TMath::Sqrt(par1.GetCovariance(0,0)) << " " << TMath::Sqrt(par1.GetCovariance(1,1)) << endl;
1082 
1083  Int_t iabs = 0, ista1 = 1;
1084  for (it1 = fVectors[ista1].begin(); it1 != fVectors[ista1].end(); ++it1) {
1085  CbmMuchTrack* tr2 = it1->second;
1086  FairTrackParam par2 = *tr2->GetParamFirst();
1087  //cout << " MUCH " << endl;
1088  //par2.Print();
1089  //cout << TMath::Sqrt(par2.GetCovariance(0,0)) << " " << TMath::Sqrt(par2.GetCovariance(1,1)) << endl;
1090  TMatrixFSym w2(5);
1091  par2.CovMatrix(w2);
1092  //w2.Invert(); // weight matrix
1093  TMatrixFSym w20 = w2;
1094  Float_t pars2[5] = {(Float_t) par2.GetX(),
1095  (Float_t) par2.GetY(),
1096  (Float_t) par2.GetTx(),
1097  (Float_t) par2.GetTy(),
1098  1.0};
1099  TMatrixF p2(5, 1, pars2);
1100  TMatrixF wp2(w2, TMatrixF::kTransposeMult, p2);
1101  wp2 += wp1;
1102  w2 += w1;
1103  w2.Invert();
1104  TMatrixF p12(w2, TMatrixF::kTransposeMult, wp2);
1105  //p12.Print();
1106 
1107  // Compute Chi2
1108  TMatrixF p122 = p12;
1109  p12 -= p1;
1110  TMatrixF wDp1(w1, TMatrixF::kMult, p12);
1111  TMatrixF chi21(p12, TMatrixF::kTransposeMult, wDp1);
1112  p122 -= p2;
1113  TMatrixF wDp2(w20, TMatrixF::kMult, p122);
1114  TMatrixF chi22(p122, TMatrixF::kTransposeMult, wDp2);
1115  Double_t c2 = chi21(0, 0) + chi22(0, 0);
1116  //cout << " Chi2: " << chi21(0,0) << " " << chi22(0,0) << " " << c2 << endl;
1117  if (c2 < 0 || c2 > fCutChi2[iabs]) continue;
1118  // Add track
1119  if (c2 < c2min) c2min = c2;
1120  //AddTrack(0, tr1, tr2, it->first, it1->first, parOk, c2); // add track
1121  } // for (it1 = fVectors[ista1].begin();
1122  if (c2min / 4 > 5) {
1123  gLogger->Info(
1124  MESSAGE_ORIGIN, "Stat.1: removed track: c2min %f", c2min / 4);
1125  fTrackArray->RemoveAt(-it->first - 1);
1126  } else
1127  cout << " Chi2: " << c2min / 4 << endl;
1128  } // for (it = fVectors[0].begin();
1129  fTrackArray->Compress();
1130 }
1131 // -------------------------------------------------------------------------
1132 
1133 // ----- Private method AddTrack1 --------------------------------------
1134 void CbmMuchMergeVectors::AddTrack1(Int_t ista0,
1135  CbmMuchTrack* tr1,
1136  CbmMuchTrack* tr2,
1137  Int_t indx1,
1138  Int_t indx2,
1139  FairTrackParam& parOk,
1140  Double_t c2) {
1141  // Store merged vector (CbmMuchTracks) into TClonesArray
1142 
1143  Int_t ntrs = fTrackArray->GetEntriesFast();
1144 
1145  CbmMuchTrack* track = new ((*fTrackArray)[ntrs]) CbmMuchTrack();
1146  track->SetParamFirst(&parOk);
1147  track->SetChiSq(c2 + tr1->GetChiSq() + tr2->GetChiSq());
1148  track->SetNDF(4 + tr1->GetNDF() + tr2->GetNDF());
1149  track->SetUniqueID(tr1->GetUniqueID()); // station number
1150  //track->SetPreviousTrackId(tr1->GetPreviousTrackId());
1151  track->SetPreviousTrackId(indx2);
1152  // Add vectors
1153  Int_t nmerged = tr1->GetNofHits();
1154  for (Int_t j = 0; j < nmerged; ++j) {
1155  if (j == 0)
1156  track->AddHit(tr1->GetHitIndex(j), tr1->GetHitType(j)); // STS
1157  else if (j == 1)
1158  track->AddHit(indx2, kMUCHSTRAWHIT); // add stat.1 vector index
1159  else
1160  track->AddHit(tr1->GetHitIndex(j - 1), tr1->GetHitType(j - 1));
1161  }
1162  if (tr1->GetFlag() == tr2->GetFlag())
1163  track->SetFlag(tr1->GetFlag());
1164  else
1165  track->SetFlag(-1);
1166 
1167  Info("AddTrack",
1168  "trID1=%i, trID2=%i, chi2=%f, merged vectors %i",
1169  tr1->GetFlag(),
1170  tr2->GetFlag(),
1171  track->GetChiSq(),
1172  track->GetNofHits());
1173 }
1174 // -------------------------------------------------------------------------
1175 
1176 ClassImp(CbmMuchMergeVectors);
CbmMuchMergeVectors.h
CbmMuchModule.h
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmTrack::GetChiSq
Double_t GetChiSq() const
Definition: CbmTrack.h:58
CbmTrack::GetParamLast
const FairTrackParam * GetParamLast() const
Definition: CbmTrack.h:62
CbmMuchMergeVectors::fRmin
Double_t fRmin[9]
Definition: CbmMuchMergeVectors.h:83
CbmMuchGeoScheme
Definition: CbmMuchGeoScheme.h:43
CbmMuchGeoScheme::GetModule
CbmMuchModule * GetModule(Int_t iStation, Int_t iLayer, Bool_t iSide, Int_t iModule) const
Definition: CbmMuchGeoScheme.cxx:238
CbmMuchStation::GetLayer
CbmMuchLayer * GetLayer(Int_t iLayer) const
Definition: CbmMuchStation.h:58
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmMuchMergeVectors::fSina
Double_t fSina[fgkPlanes]
Definition: CbmMuchMergeVectors.h:80
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmTrack.h:53
CbmMuchMergeVectors::fPoints
TClonesArray * fPoints
Definition: CbmMuchMergeVectors.h:66
CbmTrack::GetFlag
Int_t GetFlag() const
Definition: CbmTrack.h:57
CbmMuchMergeVectors::fVectors
std::map< Int_t, CbmMuchTrack * > fVectors[fgkStat]
Definition: CbmMuchMergeVectors.h:77
CbmMuchMergeVectors::fZabs0
Double_t fZabs0[9][2]
Definition: CbmMuchMergeVectors.h:87
lun1
FILE * lun1
Definition: CbmMuchMergeVectors.cxx:39
CbmMuchMergeVectors::MatchVectors
void MatchVectors()
Definition: CbmMuchMergeVectors.cxx:328
CbmMuchMergeVectors::fMatr
std::map< Int_t, TMatrixDSym * > fMatr
Definition: CbmMuchMergeVectors.h:86
CbmMuchMergeVectors::fHits
TClonesArray * fHits
Definition: CbmMuchMergeVectors.h:64
CbmKFTrack::GetTrackParam
void GetTrackParam(FairTrackParam &track)
Definition: CbmKFTrack.cxx:45
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMuchMergeVectors::PassAbsorber
void PassAbsorber(Int_t ist, Double_t *zabs, Double_t x0, FairTrackParam &parFirst, TMatrixFSym &cov, Int_t pFlag=0)
Definition: CbmMuchMergeVectors.cxx:658
CbmMuchMergeVectors::fNdoubl
Int_t fNdoubl
Definition: CbmMuchMergeVectors.h:74
CbmMuchMergeVectors::fRmax
Double_t fRmax[9]
Definition: CbmMuchMergeVectors.h:84
CbmMuchMergeVectors::AddTrack
void AddTrack(Int_t ista0, CbmMuchTrack *tr1, CbmMuchTrack *tr2, Int_t indx1, Int_t indx2, FairTrackParam &parOk, Double_t c2)
Definition: CbmMuchMergeVectors.cxx:741
CbmTrack::SetPreviousTrackId
void SetPreviousTrackId(Int_t previousTrackId)
Definition: CbmTrack.h:72
CbmMuchMergeVectors::CbmMuchMergeVectors
CbmMuchMergeVectors()
Definition: CbmMuchMergeVectors.cxx:42
CbmTrack::SetParamLast
void SetParamLast(const FairTrackParam *par)
Definition: CbmTrack.h:76
CbmTrack::SetNDF
void SetNDF(Int_t ndf)
Definition: CbmTrack.h:71
CbmMuchLayer::GetSide
CbmMuchLayerSide * GetSide(Bool_t side)
Definition: CbmMuchLayer.h:49
CbmMuchTrack
Definition: CbmMuchTrack.h:16
CbmHit::GetPlaneId
virtual Int_t GetPlaneId() const
Definition: CbmHit.h:97
CbmStsTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmStsTrack.h:76
CbmTrack::SetFlag
void SetFlag(Int_t flag)
Definition: CbmTrack.h:69
CbmMuchMergeVectors::~CbmMuchMergeVectors
virtual ~CbmMuchMergeVectors()
Definition: CbmMuchMergeVectors.cxx:66
CbmMuchMergeVectors::fTrStsMatch
TClonesArray * fTrStsMatch
Definition: CbmMuchMergeVectors.h:70
CbmMuchGeoScheme::GetNStations
Int_t GetNStations() const
Definition: CbmMuchGeoScheme.h:79
CbmStsTrack.h
Data class for STS tracks.
CbmMuchMergeVectors::Init
virtual InitStatus Init()
Definition: CbmMuchMergeVectors.cxx:76
CbmMuchPoint.h
CbmTrack::GetMatch
CbmMatch * GetMatch() const
Definition: CbmTrack.h:63
CbmMuchTrack.h
CbmMuchMergeVectors::fgkPlanes
static const Int_t fgkPlanes
Definition: CbmMuchMergeVectors.h:56
CbmMuchMergeVectors::fVecArray
TClonesArray * fVecArray
Definition: CbmMuchMergeVectors.h:68
CbmMuchMergeVectors::fCosa
Double_t fCosa[fgkPlanes]
Definition: CbmMuchMergeVectors.h:79
CbmMuchMergeVectors
Definition: CbmMuchMergeVectors.h:23
CbmMuchStation::GetRmax
Double_t GetRmax() const
Definition: CbmMuchStation.h:52
CbmMuchMergeVectors::Exec
virtual void Exec(Option_t *opt)
Definition: CbmMuchMergeVectors.cxx:157
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmTrackMatchNew.h
ClassImp
ClassImp(CbmMuchMergeVectors)
CbmTrack::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTrack.h:70
CbmMuchMergeVectors::fCutChi2
Double_t fCutChi2[9]
Definition: CbmMuchMergeVectors.h:81
CbmMuchStation::GetDetectorId
Int_t GetDetectorId() const
Definition: CbmMuchStation.h:48
CbmTrack::GetNDF
Int_t GetNDF() const
Definition: CbmTrack.h:59
CbmMuchLayerSide
Definition: CbmMuchLayerSide.h:22
kHIT
@ kHIT
Definition: CbmHit.h:17
CbmKFTrack.h
CbmMuchMergeVectors::fStatFirst
Int_t fStatFirst
Definition: CbmMuchMergeVectors.h:75
CbmTrack::GetPreviousTrackId
Int_t GetPreviousTrackId() const
Definition: CbmTrack.h:60
CbmMuchMergeVectors::fTracksSts
TClonesArray * fTracksSts
Definition: CbmMuchMergeVectors.h:69
CbmTrack::GetHitType
HitType GetHitType(Int_t iHit) const
Definition: CbmTrack.h:55
CbmMuchMergeVectors::GetVectors
void GetVectors()
Definition: CbmMuchMergeVectors.cxx:196
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmMuchModule
Definition: CbmMuchModule.h:24
CbmMuchFindVectorsGem::GetAbsorbers
Int_t GetAbsorbers(Double_t zabs[9][2], Double_t *x0abs)
Definition: CbmMuchFindVectorsGem.cxx:1356
CbmMuchAddress::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchAddress.h:103
CbmMuchMergeVectors::fX0abs
Double_t fX0abs[9]
Definition: CbmMuchMergeVectors.h:88
CbmMuchMergeVectors::fTrackArray
TClonesArray * fTrackArray
Definition: CbmMuchMergeVectors.h:62
CbmMuchMergeVectors::Finish
virtual void Finish()
Definition: CbmMuchMergeVectors.cxx:186
CbmMuchModule::GetDetectorType
Int_t GetDetectorType() const
Definition: CbmMuchModule.h:52
CbmMuchMergeVectors::SetParContainers
virtual void SetParContainers()
Definition: CbmMuchMergeVectors.cxx:153
CbmMuchMergeVectors::fZ0
Double_t fZ0[9]
Definition: CbmMuchMergeVectors.h:82
CbmMuchFindVectorsGem.h
CbmHit
Definition: CbmHit.h:38
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmMuchMergeVectors::fTracksLit
TClonesArray * fTracksLit
Definition: CbmMuchMergeVectors.h:71
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
kMUCHSTRAWHIT
@ kMUCHSTRAWHIT
Definition: CbmHit.h:24
CbmTrack::SetParamFirst
void SetParamFirst(const FairTrackParam *par)
Definition: CbmTrack.h:75
hits
static vector< vector< QAHit > > hits
Definition: CbmTofHitFinderTBQA.cxx:114
CbmMuchMergeVectors::AddTrack1
void AddTrack1(Int_t ista0, CbmMuchTrack *tr1, CbmMuchTrack *tr2, Int_t indx1, Int_t indx2, FairTrackParam &parOk, Double_t c2)
Definition: CbmMuchMergeVectors.cxx:1134
CbmMuchGeoScheme::GetStation
CbmMuchStation * GetStation(Int_t iStation) const
Definition: CbmMuchGeoScheme.cxx:209
CbmMuchStation::GetNLayers
Int_t GetNLayers() const
Definition: CbmMuchStation.h:50
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmMuchDigiMatch.h
CbmMuchMergeVectors::AddStation1
void AddStation1()
Definition: CbmMuchMergeVectors.cxx:1015
CbmMuchMergeVectors::fDz
Double_t fDz[fgkPlanes]
Definition: CbmMuchMergeVectors.h:78
CbmMuchMergeVectors::MergeVectors
void MergeVectors()
Definition: CbmMuchMergeVectors.cxx:473
CbmMuchMergeVectors::fGemHits
TClonesArray * fGemHits
Definition: CbmMuchMergeVectors.h:65
CbmMuchLayer
Definition: CbmMuchLayer.h:21
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmKFTrack
Definition: CbmKFTrack.h:21
CbmMuchLayerSide::GetZ
Double_t GetZ()
Definition: CbmMuchLayerSide.h:49
CbmMuchStation::GetRmin
Double_t GetRmin() const
Definition: CbmMuchStation.h:51
CbmMuchStation.h
CbmTrack::AddHit
void AddHit(Int_t index, HitType type)
Definition: CbmTrack.cxx:75
CbmMuchFindVectorsGem
Definition: CbmMuchFindVectorsGem.h:27
CbmMuchMergeVectors::RemoveClones
void RemoveClones(Int_t ibeg, Int_t iabs, std::map< Int_t, CbmMuchTrack * > &mapMerged)
Definition: CbmMuchMergeVectors.cxx:808
CbmMuchMergeVectors::fgkStat
static const Int_t fgkStat
Definition: CbmMuchMergeVectors.h:57
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0)
Access to i-th hit.
Definition: CbmKFTrackInterface.cxx:39
CbmMuchMergeVectors::SelectTracks
void SelectTracks()
Definition: CbmMuchMergeVectors.cxx:895
CbmMuchMergeVectors::fGeoScheme
CbmMuchGeoScheme * fGeoScheme
Definition: CbmMuchMergeVectors.h:61