CbmRoot
riplet/LxTrackAna.cxx
Go to the documentation of this file.
1 #include "LxTrackAna.h"
2 #include "CbmKFParticle.h"
3 #include "CbmKFTrack.h"
4 #include "CbmMCTrack.h"
5 #include "CbmMuchCluster.h"
6 #include "CbmMuchDigiMatch.h"
7 #include "CbmMuchGeoScheme.h"
8 #include "CbmMuchPoint.h"
9 #include "CbmStsAddress.h"
10 #include "CbmStsPoint.h"
11 #include "CbmStsTrack.h"
12 #include "TDatabasePDG.h"
13 #include "TH1.h"
14 #include "TH2.h"
15 #include <dirent.h>
16 #include <iostream>
17 #include <sys/stat.h>
18 #include <sys/types.h>
19 
20 ClassImp(LxTrackAnaTriplet)
21 
22  using namespace std;
23 
24 //static Double_t xRms = 1.005;// Averaged MC
25 static Double_t xRms = 1.202; // Nearest hit
26 static Double_t xRms2 = xRms * xRms;
27 //static Double_t yRms = 0.9467;// Averaged MC
28 static Double_t yRms = 1.061; // Nearest hit
29 static Double_t yRms2 = yRms * yRms;
30 //static Double_t txRms = 0.02727;// Averaged MC
31 static Double_t txRms = 0.02426; // Nearest hit
32 static Double_t txRms2 = txRms * txRms;
33 //static Double_t tyRms = 0.01116;// Averaged MC
34 static Double_t tyRms = 0.01082; // Nearest hit
35 static Double_t tyRms2 = tyRms * tyRms;
36 static Double_t cutCoeff = 3.0;
37 
38 static TH1F* muchStsBreakX = 0;
39 static TH1F* muchStsBreakY = 0;
40 static TH1F* muchStsBreakTx = 0;
41 static TH1F* muchStsBreakTy = 0;
42 
43 static TH1F* stsMuchBreakX = 0;
44 static TH1F* stsMuchBreakY = 0;
45 
46 static TH1F* signalChi2 = 0;
47 static TH1F* bgrChi2 = 0;
48 
49 static TH1F* bgrInvMass = 0;
50 static list<LxSimpleTrack*> positiveTracks;
51 static list<LxSimpleTrack*> negativeTracks;
52 static TH1F* sigInvMass = 0;
53 
54 static TH1F* nearestHitDist[LXSTATIONS] = {0};
55 static TH1F* hitsDist[LXSTATIONS] = {0};
56 
57 static TH1F* muPlusStsTxDiff = 0;
58 static TH1F* muMinusStsTxDiff = 0;
59 static TH1F* muPlusStsXDiff = 0;
60 static TH1F* muMinusStsXDiff = 0;
61 static TH1F* muPlusVertexTxDiff = 0;
62 static TH1F* muMinusVertexTxDiff = 0;
63 
64 static TH2F* muPlusStsBeginTxDiff2D = 0;
65 static TH2F* muMinusStsBeginTxDiff2D = 0;
66 
67 static TH1F* deltaPhiPi = 0;
68 
69 static TH1F* jPsiMuonsMomsHisto = 0;
70 
71 struct MomVsTxRange {
72  Double_t momLow;
73  Double_t momHigh;
74  Double_t txLow;
75  Double_t txHigh;
76 };
77 
78 static bool momFitTxBreak(Double_t mom, Double_t txBreak) {
79  if (mom < 3.0) return false;
80 
81  if (txBreak < 0) txBreak = -txBreak;
82 
83  Double_t inv = mom * txBreak;
84  return inv > 0.18 && inv < 0.52;
85 }
86 
87 void LxSimpleTrack::RebindMuchTrack() {
88  linkedStsTrack = 0;
89 
90  while (!linkedStsTracks.empty()) {
91  pair<LxSimpleTrack*, Double_t> trackDesc = linkedStsTracks.front();
92  linkedStsTracks.pop_front();
93  LxSimpleTrack* anotherMuchTrack = trackDesc.first->linkedMuchTrack.first;
94 
95  if (0 == anotherMuchTrack
96  || trackDesc.second < trackDesc.first->linkedMuchTrack.second) {
97  trackDesc.first->linkedMuchTrack.first = this;
98  trackDesc.first->linkedMuchTrack.second = trackDesc.second;
99  linkedStsTrack = trackDesc.first;
100 
101  if (0 != anotherMuchTrack) anotherMuchTrack->RebindMuchTrack();
102 
103  break;
104  }
105  }
106 }
107 
108 LxTrackAnaTriplet::LxTrackAnaTriplet()
109  : listMCTracks(0)
110  , listStsPts(0)
111  , listMuchPts(0)
112  , listMuchPixelHits(0)
113  , listMuchClusters(0)
114  , listMuchPixelDigiMatches(0)
115  , allTracks()
116  , posTracks()
117  , negTracks()
118  , superEventTracks(0)
119  , superEventBrachTrack(0, 0, 0, 0, 0, 0, 0, 0)
120  , useHitsInStat(false)
121  , averagePoints(false)
122  , dontTouchNonPrimary(true)
123  , useChargeSignInCuts(false)
124  , buildConnectStat(false)
125  , buildBgrInvMass(false)
126  , buildSigInvMass(false)
127  , joinData(false)
128  , buildNearestHitDist(false)
129  , cropHits(false)
130  , buildSegmentsStat(true)
131  , particleType("jpsi")
132  , segmentsAnalyzer(*this) {}
133 
134 LxTrackAnaTriplet::~LxTrackAnaTriplet() { Clean(); }
135 
136 void LxTrackAnaTriplet::Clean() {
137  for (vector<LxSimpleTrack*>::iterator i = allTracks.begin();
138  i != allTracks.end();
139  ++i)
140  delete *i;
141 
142  allTracks.clear();
143  posTracks.clear();
144  negTracks.clear();
145 }
146 
147 InitStatus LxTrackAnaTriplet::Init() {
148  FairRootManager* fManager = FairRootManager::Instance();
149  listMCTracks = static_cast<TClonesArray*>(fManager->GetObject("MCTrack"));
150 
151  if (0 == listMCTracks) return kFATAL;
152 
153  listStsPts = static_cast<TClonesArray*>(fManager->GetObject("StsPoint"));
154 
155  if (0 == listStsPts) return kFATAL;
156 
157  listMuchPts = static_cast<TClonesArray*>(fManager->GetObject("MuchPoint"));
158 
159  if (0 == listMuchPts) return kFATAL;
160 
161  if (useHitsInStat) {
162  listMuchPixelHits =
163  static_cast<TClonesArray*>(fManager->GetObject("MuchPixelHit"));
164 
165  if (0 == listMuchPixelHits) return kFATAL;
166 
167  listMuchClusters =
168  static_cast<TClonesArray*>(fManager->GetObject("MuchCluster"));
169 
170  if (0 == listMuchClusters) return kFATAL;
171 
172  listMuchPixelDigiMatches =
173  static_cast<TClonesArray*>(fManager->GetObject("MuchDigiMatch"));
174 
175  if (0 == listMuchPixelDigiMatches) return kFATAL;
176  }
177 
178  if (buildConnectStat) {
179  muchStsBreakX = new TH1F(
180  "muchStsBreakX", "Break in prediction of X in STS", 100, -20., 20.);
181  muchStsBreakX->StatOverflows();
182  muchStsBreakY = new TH1F(
183  "muchStsBreakY", "Break in prediction of Y in STS", 100, -20., 20.);
184  muchStsBreakY->StatOverflows();
185  muchStsBreakTx = new TH1F(
186  "muchStsBreakTx", "Break in prediction of Tx in STS", 100, -0.15, 0.15);
187  muchStsBreakTx->StatOverflows();
188  muchStsBreakTy = new TH1F(
189  "muchStsBreakTy", "Break in prediction of Ty in STS", 100, -0.15, 0.15);
190  muchStsBreakTy->StatOverflows();
191 
192  stsMuchBreakX = new TH1F(
193  "stsMuchBreakX", "Break in prediction of X in MUCH", 100, -20., 20.);
194  stsMuchBreakX->StatOverflows();
195  stsMuchBreakY = new TH1F(
196  "stsMuchBreakY", "Break in prediction of Y in MUCH", 100, -20., 20.);
197  stsMuchBreakY->StatOverflows();
198 
199  signalChi2 = new TH1F("signalChi2", "Chi2 of signal", 100, 0., 15.);
200  signalChi2->StatOverflows();
201  bgrChi2 = new TH1F("bgrChi2", "Chi2 of background", 100, 0., 20.);
202  bgrChi2->StatOverflows();
203  } // if (buildConnectStat)
204 
205  if (buildBgrInvMass) {
206  if (joinData) {
207  bgrInvMass = new TH1F("bgrInvMass",
208  "Invariant mass distribution for background",
209  1000,
210  2.,
211  4.);
212  bgrInvMass->StatOverflows();
213  } else {
214  superEventTracks =
215  new TTree("SuperEventTracks", "Tracks for building a super event");
216  superEventTracks->Branch(
217  "tracks", &superEventBrachTrack.px, "px/D:py:pz:e:charge");
218  }
219  } // if (buildBgrInvMass)
220 
221  if (buildSigInvMass) {
222  sigInvMass = new TH1F(
223  "sigInvMass", "Invariant mass distribution for signal", 1000, 2., 4.);
224  sigInvMass->StatOverflows();
225  } // if (buildSigInvMass)
226 
227  if (buildNearestHitDist && useHitsInStat) {
228  char name[32];
229  char title[128];
230 
231  for (Int_t i = 0; i < LXSTATIONS; ++i) {
232  sprintf(name, "nearestHitDist_%d", i);
233  sprintf(
234  title, "Distance from a MC point to the nearest hit at %d station", i);
235  nearestHitDist[i] = new TH1F(name, title, 100, 0., 5.);
236  nearestHitDist[i]->StatOverflows();
237  sprintf(name, "hitsDist_%d", i);
238  sprintf(title, "Distance from a MC point to the hits at %d station", i);
239  hitsDist[i] = new TH1F(name, title, 100, 0., 7.);
240  hitsDist[i]->StatOverflows();
241  }
242  }
243 
244  muPlusStsTxDiff =
245  new TH1F("muPlusStsTxDiff", "muPlusStsTxDiff", 100, -0.2, 0.2);
246  muPlusStsTxDiff->StatOverflows();
247  muMinusStsTxDiff =
248  new TH1F("muMinusStsTxDiff", "muMinusStsTxDiff", 100, -0.2, 0.2);
249  muMinusStsTxDiff->StatOverflows();
250  muPlusStsXDiff =
251  new TH1F("muPlusStsXDiff", "muPlusStsXDiff", 100, -10.0, 10.0);
252  muPlusStsXDiff->StatOverflows();
253  muMinusStsXDiff =
254  new TH1F("muMinusStsXDiff", "muMinusStsXDiff", 100, -10.0, 10.0);
255  muMinusStsXDiff->StatOverflows();
256  muPlusVertexTxDiff =
257  new TH1F("muPlusVertexTxDiff", "muPlusVertexTxDiff", 100, -0.2, 0.2);
258  muPlusVertexTxDiff->StatOverflows();
259  muMinusVertexTxDiff =
260  new TH1F("muMinusVertexTxDiff", "muMinusVertexTxDiff", 100, -0.2, 0.2);
261  muMinusVertexTxDiff->StatOverflows();
262 
263  muPlusStsBeginTxDiff2D = new TH2F("muPlusStsBeginTxDiff2D",
264  "muPlusStsBeginTxDiff2D",
265  100,
266  0.,
267  25.,
268  100,
269  -0.2,
270  0.2);
271  muPlusStsBeginTxDiff2D->StatOverflows();
272  muMinusStsBeginTxDiff2D = new TH2F("muMinusStsBeginTxDiff2D",
273  "muMinusStsBeginTxDiff2D",
274  100,
275  0.,
276  25.0,
277  100,
278  -0.2,
279  0.2);
280  muMinusStsBeginTxDiff2D->StatOverflows();
281 
282  deltaPhiPi = new TH1F("deltaPhiPi", "deltaPhiPi", 100, 0., 1.0);
283  deltaPhiPi->StatOverflows();
284 
285  jPsiMuonsMomsHisto = new TH1F(
286  "jPsiMuonsMomsHisto", "J/Psi muons momenta distribution", 200, 0., 25.);
287  jPsiMuonsMomsHisto->StatOverflows();
288 
289  segmentsAnalyzer.Init();
290 
291  return kSUCCESS;
292 }
293 
294 static void SaveHisto(TH1* histo, const char* particleType, const char* name) {
295  if (!saveHistos) return;
296 
297  char dir_name[256];
298  sprintf(dir_name, "configuration.%s", particleType);
299  DIR* dir = opendir(dir_name);
300 
301  if (dir)
302  closedir(dir);
303  else
304  mkdir(dir_name, 0700);
305 
306  char file_name[256];
307  sprintf(file_name, "%s/%s", dir_name, name);
308  TFile fh(file_name, "RECREATE");
309  histo->Write();
310  fh.Close();
311  delete histo;
312 }
313 
314 static void BuildInvMass(list<LxSimpleTrack*>& pTracks,
315  list<LxSimpleTrack*>& nTracks,
316  TH1* histo) {
317  for (list<LxSimpleTrack*>::iterator i = pTracks.begin(); i != pTracks.end();
318  ++i) {
319  LxSimpleTrack* posTrack = *i;
320 
321  for (list<LxSimpleTrack*>::iterator j = nTracks.begin(); j != nTracks.end();
322  ++j) {
323  LxSimpleTrack* negTrack = *j;
324  Double_t E12 = posTrack->e + negTrack->e;
325  Double_t px12 = posTrack->px + negTrack->px;
326  Double_t py12 = posTrack->py + negTrack->py;
327  Double_t pz12 = posTrack->pz + negTrack->pz;
328  Double_t m122 = E12 * E12 - px12 * px12 - py12 * py12 - pz12 * pz12;
329  Double_t m12 = m122 > 1.e-20 ? sqrt(m122) : 0.;
330  histo->Fill(m12);
331  }
332  }
333 }
334 
335 static void BuildInvMass2(list<CbmStsTrack*>& stsTracks, TH1* /*histo*/) {
336  for (list<CbmStsTrack*>::iterator i = stsTracks.begin(); i != stsTracks.end();
337  ++i) {
338  CbmStsTrack* posTrack = *i;
339  //CbmStsTrack t1(*posTrack);
340  //extFitter.DoFit(&t1, 13);
341  //Double_t chi2Prim = extFitter.GetChiToVertex(&t1, fPrimVtx);
342  const FairTrackParam* t1param = posTrack->GetParamFirst();
343  //extFitter.Extrapolate(&t1, fPrimVtx->GetZ(), &t1param);
344  CbmKFTrack muPlus(*posTrack);
345  Double_t t1Qp = t1param->GetQp();
346  list<CbmStsTrack*>::iterator j = i;
347  ++j;
348 
349  for (; j != stsTracks.end(); ++j) {
350  CbmStsTrack* negTrack = *j;
351  //CbmStsTrack t2(*negTrack);
352  //extFitter.DoFit(&t2, 13);
353  //chi2Prim = extFitter.GetChiToVertex(&t2, fPrimVtx);
354  const FairTrackParam* t2param = negTrack->GetParamLast();
355  //extFitter.Extrapolate(&t2, fPrimVtx->GetZ(), &t2param);
356  Double_t t2Qp = t2param->GetQp();
357 
358  if (t1Qp * t2Qp >= 0) continue;
359 
360  CbmKFTrack muMinus(*negTrack);
361  vector<CbmKFTrackInterface*> kfData;
362 
363  if (t1Qp > 0) {
364  kfData.push_back(&muPlus);
365  kfData.push_back(&muMinus);
366  } else {
367  kfData.push_back(&muMinus);
368  kfData.push_back(&muPlus);
369  }
370 
371  /*CbmKFParticle DiMu;
372  DiMu.Construct(kfData, 0);
373  DiMu.TransportToDecayVertex();
374  Double_t m, merr;
375  DiMu.GetMass(m, merr);
376  histo->Fill(m);*/
377  }
378  }
379 }
380 
381 void LxTrackAnaTriplet::FinishTask() {
382  TFile* curFile = TFile::CurrentFile();
383 
384  if (buildConnectStat) {
385  SaveHisto(muchStsBreakX, particleType.Data(), "muchStsBreakX.root");
386  SaveHisto(muchStsBreakY, particleType.Data(), "muchStsBreakY.root");
387  SaveHisto(muchStsBreakTx, particleType.Data(), "muchStsBreakTx.root");
388  SaveHisto(muchStsBreakTy, particleType.Data(), "muchStsBreakTy.root");
389 
390  SaveHisto(stsMuchBreakX, particleType.Data(), "stsMuchBreakX.root");
391  SaveHisto(stsMuchBreakY, particleType.Data(), "stsMuchBreakY.root");
392 
393  SaveHisto(signalChi2, particleType.Data(), "signalChi2.root");
394  SaveHisto(bgrChi2, particleType.Data(), "bgrChi2.root");
395  } // if (buildConnectStat)
396 
397  if (buildBgrInvMass) {
398  if (joinData) {
399  TFile fh("tracks_tree.root");
400  superEventTracks = static_cast<TTree*>(fh.Get("SuperEventTracks"));
401  //LxSimpleTrack st(0, 0, 0, 0, 0, 0, 0, 0);
402  //superEventTracks->SetBranchAddress("tracks", &st.px);
403  CbmStsTrack* st = new CbmStsTrack;
404  superEventTracks->SetBranchAddress("tracks", &st);
405  list<CbmStsTrack*> stsTracks;
406  Int_t nEnt = superEventTracks->GetEntries();
407 
408  for (Int_t i = 0; i < nEnt; ++i) {
409  superEventTracks->GetEntry(i);
410  //LxSimpleTrack* t = new LxSimpleTrack(0, 0, 0, 0, st.px, st.py, st.pz, st.e);
411  //t->charge = st.charge;
412  CbmStsTrack* t = new CbmStsTrack(*st);
413  stsTracks.push_back(t);
414  }
415 
416  BuildInvMass2(stsTracks, bgrInvMass);
417  SaveHisto(bgrInvMass, particleType.Data(), "bgrInvMass.root");
418 
419  for (list<CbmStsTrack*>::iterator i = stsTracks.begin();
420  i != stsTracks.end();
421  ++i)
422  delete *i;
423  } else {
424  TFile fh("tracks_tree.root", "RECREATE");
425  superEventTracks->Write();
426  fh.Close();
427  delete superEventTracks;
428  }
429  } // if (buildBgrInvMass)
430 
431  if (buildSigInvMass)
432  SaveHisto(sigInvMass, particleType.Data(), "sigInvMass.root");
433 
434  if (buildNearestHitDist && useHitsInStat) {
435  char fileName[32];
436 
437  for (Int_t i = 0; i < LXSTATIONS; ++i) {
438  sprintf(fileName, "%s.root", nearestHitDist[i]->GetName());
439  SaveHisto(nearestHitDist[i], particleType.Data(), fileName);
440  sprintf(fileName, "%s.root", hitsDist[i]->GetName());
441  SaveHisto(hitsDist[i], particleType.Data(), fileName);
442  }
443  }
444 
445  SaveHisto(muPlusStsTxDiff, particleType.Data(), "muPlusStsTxDiff.root");
446  SaveHisto(muMinusStsTxDiff, particleType.Data(), "muMinusStsTxDiff.root");
447  SaveHisto(muPlusStsXDiff, particleType.Data(), "muPlusStsXDiff.root");
448  SaveHisto(muMinusStsXDiff, particleType.Data(), "muMinusStsXDiff.root");
449  SaveHisto(muPlusVertexTxDiff, particleType.Data(), "muPlusVertexTxDiff.root");
450  SaveHisto(
451  muMinusVertexTxDiff, particleType.Data(), "muMinusVertexTxDiff.root");
452 
453  SaveHisto(
454  muPlusStsBeginTxDiff2D, particleType.Data(), "muPlusStsBeginTxDiff2D.root");
455  SaveHisto(muMinusStsBeginTxDiff2D,
456  particleType.Data(),
457  "muMinusStsBeginTxDiff2D.root");
458 
459  SaveHisto(deltaPhiPi, particleType.Data(), "deltaPhiPi.root");
460 
461  SaveHisto(jPsiMuonsMomsHisto, particleType.Data(), "jPsiMuonsMomsHisto.root");
462 
463  segmentsAnalyzer.Finish();
464 
465  TFile::CurrentFile() = curFile;
466  FairTask::FinishTask();
467 }
468 
469 // Our goal here is to investigate various properties of Monte Carlo tracks derivable from points of there intersections
470 // with detector stations. At the same time in MUCH we use not the Monte Carlo points but hits corresponding to them.
471 // -- This is done to make the statistical properties more realistic.
472 void LxTrackAnaTriplet::Exec(Option_t*) {
473  Clean();
474 
475  Int_t nEnt = listMCTracks->GetEntries();
476  cout << "There are: " << nEnt << " of MC tracks" << endl;
477 
478  for (Int_t i = 0; i < nEnt; ++i) {
479  CbmMCTrack* mct = static_cast<CbmMCTrack*>(listMCTracks->At(i));
480  LxSimpleTrack* t = new LxSimpleTrack(mct->GetPdgCode(),
481  mct->GetMotherId(),
482  mct->GetP(),
483  mct->GetPt(),
484  mct->GetPx(),
485  mct->GetPy(),
486  mct->GetPz(),
487  mct->GetEnergy());
488 
489  if (t->motherId >= 0) t->parent = allTracks[t->motherId];
490 
491  allTracks.push_back(t);
492  }
493 
494  nEnt = listStsPts->GetEntries();
495  cout << "There are: " << nEnt << " of STS MC points" << endl;
496 
497  for (Int_t i = 0; i < nEnt; ++i) {
498  CbmStsPoint* stsPt = static_cast<CbmStsPoint*>(listStsPts->At(i));
499  Int_t mcTrackId = stsPt->GetTrackID();
500  LxSimpleTrack* track = allTracks[mcTrackId];
501  TVector3 xyzI, xyzO;
502  stsPt->Position(xyzI);
503  stsPt->PositionOut(xyzO);
504  TVector3 xyz = .5 * (xyzI + xyzO);
505  TVector3 pxyzI, pxyzO;
506  stsPt->Momentum(pxyzI);
507  stsPt->MomentumOut(pxyzO);
508  TVector3 pxyz = .5 * (pxyzI + pxyzO);
509  LxSimplePoint point(
510  xyz.X(), xyz.Y(), xyz.Z(), pxyz.X() / pxyz.Z(), pxyz.Y() / pxyz.Z());
511  Int_t stationNr =
512  CbmStsAddress::GetElementId(stsPt->GetDetectorID(), kStsStation);
513  track->stsPoints[stationNr].push_back(point);
514  }
515 
516  nEnt = listMuchPts->GetEntries();
517  cout << "There are: " << nEnt << " of MUCH MC points" << endl;
518 
519  for (Int_t i = 0; i < nEnt; ++i) {
520  CbmMuchPoint* mcPoint = static_cast<CbmMuchPoint*>(listMuchPts->At(i));
521  Int_t mcTrackId = mcPoint->GetTrackID();
522  LxSimpleTrack* track = allTracks[mcTrackId];
523  Int_t stationNr =
524  CbmMuchGeoScheme::GetStationIndex(mcPoint->GetDetectorId());
525  Int_t layerNr = CbmMuchGeoScheme::GetLayerIndex(mcPoint->GetDetectorId());
526  TVector3 xyzI, xyzO;
527  mcPoint->Position(xyzI);
528  mcPoint->PositionOut(xyzO);
529  TVector3 xyz = .5 * (xyzI + xyzO);
530  TVector3 pxyzI, pxyzO;
531  mcPoint->Momentum(pxyzI);
532  mcPoint->MomentumOut(pxyzO);
533  TVector3 pxyz = .5 * (pxyzI + pxyzO);
534  LxSimplePoint point(
535  xyz.X(), xyz.Y(), xyz.Z(), pxyz.X() / pxyz.Z(), pxyz.Y() / pxyz.Z());
536 
537  if (useHitsInStat)
538  track->muchMCPts[stationNr][layerNr].push_back(point);
539  else
540  track->muchPoints[stationNr][layerNr].push_back(point);
541  }
542 
543  if (useHitsInStat) {
544  nEnt = listMuchPixelHits->GetEntries();
545  cout << "There are: " << nEnt << " of MUCH pixel hits" << endl;
546 
547  for (Int_t i = 0; i < nEnt; ++i) {
548  CbmMuchPixelHit* mh =
549  static_cast<CbmMuchPixelHit*>(listMuchPixelHits->At(i));
550 
551  Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(mh->GetAddress());
552  Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(mh->GetAddress());
553  TVector3 pos, err;
554  mh->Position(pos);
555  mh->PositionError(err);
556  Double_t x = pos.X();
557  Double_t y = pos.Y();
558  Double_t z = pos.Z();
559 
560  if (x != x || y != y || z != z) // Test for NaN
561  continue;
562 
563  set<LxSimpleTrack*> tracks;
564 
565  Int_t clusterId = mh->GetRefId();
566  CbmMuchCluster* cluster =
567  static_cast<CbmMuchCluster*>(listMuchClusters->At(clusterId));
568  Int_t nDigis = cluster->GetNofDigis();
569 
570  for (Int_t j = 0; j < nDigis; ++j) {
571  CbmMuchDigiMatch* digiMatch = static_cast<CbmMuchDigiMatch*>(
572  listMuchPixelDigiMatches->At(cluster->GetDigi(j)));
573  Int_t nMCs = digiMatch->GetNofLinks();
574 
575  for (Int_t k = 0; k < nMCs; ++k) {
576  const CbmLink& lnk = digiMatch->GetLink(k);
577  Int_t mcIndex = lnk.GetIndex();
578  CbmMuchPoint* mcPoint =
579  static_cast<CbmMuchPoint*>(listMuchPts->At(mcIndex));
580  Int_t mcTrackId = mcPoint->GetTrackID();
581  tracks.insert(allTracks[mcTrackId]);
582  //Int_t stationNr = CbmMuchGeoScheme::GetStationIndex(mcPoint->GetDetectorId());
583  //Int_t layerNr = CbmMuchGeoScheme::GetLayerIndex(mcPoint->GetDetectorId());
584  //LxSimplePoint point(x, y, z, 0, 0);
585  //track->muchPoints[stationNr][layerNr].push_back(point);
586  }
587  } // j
588 
589  for (set<LxSimpleTrack*>::iterator j = tracks.begin(); j != tracks.end();
590  ++j) {
591  LxSimplePoint point(x, y, z, 0, 0);
592  LxSimpleTrack* track = *j;
593  track->muchPoints[stationNumber][layerNumber].push_back(point);
594  } // j
595  } // i
596 
597  nEnt = listMuchClusters->GetEntries();
598  cout << "There are: " << nEnt << " of MUCH clusters" << endl;
599 
600  nEnt = listMuchPixelDigiMatches->GetEntries();
601  cout << "There are: " << nEnt << " of MUCH pixel digi matches" << endl;
602  } //if (useHitsInStat)
603 
604  if (averagePoints) AveragePoints();
605 
606  if (buildConnectStat) BuildStatistics();
607 
608  //Connect(false);
609  Connect(true);
610 
611  if (buildSigInvMass) BuildInvMass(posTracks, negTracks, sigInvMass);
612 
613  if (buildSegmentsStat) segmentsAnalyzer.BuildStatistics();
614 }
615 
616 static inline void AveragePoints(list<LxSimplePoint>& points) {
617  if (points.empty() || points.size() == 1) return;
618 
619  Double_t x = 0;
620  Double_t y = 0;
621  Double_t z = 0;
622  Double_t tx = 0;
623  Double_t ty = 0;
624 
625  for (list<LxSimplePoint>::iterator i = points.begin(); i != points.end();
626  ++i) {
627  LxSimplePoint p = *i;
628  x += p.x;
629  y += p.y;
630  z += p.z;
631  tx += p.tx;
632  ty += p.ty;
633  }
634 
635  x /= points.size();
636  y /= points.size();
637  z /= points.size();
638  tx /= points.size();
639  ty /= points.size();
640  LxSimplePoint averagedPoint(x, y, z, tx, ty);
641  points.clear();
642  points.push_back(averagedPoint);
643 }
644 
645 // If hits are used in statistics average MC points in MUCH.
646 static inline void AveragePoints(LxSimpleTrack* track, bool useHitsInStat) {
647  for (Int_t i = 0; i < LXSTSSTATIONS; ++i)
648  AveragePoints(track->stsPoints[i]);
649 
650  if (useHitsInStat) {
651  for (Int_t i = 0; i < LXSTATIONS; ++i) {
652  for (Int_t j = 0; j < LXLAYERS; ++j)
653  AveragePoints(track->muchMCPts[i][j]);
654  }
655  } else {
656  for (Int_t i = 0; i < LXSTATIONS; ++i) {
657  for (Int_t j = 0; j < LXLAYERS; ++j)
658  AveragePoints(track->muchPoints[i][j]);
659  }
660  }
661 }
662 
663 void LxTrackAnaTriplet::AveragePoints() {
664  for (vector<LxSimpleTrack*>::iterator i = allTracks.begin();
665  i != allTracks.end();
666  ++i)
667  ::AveragePoints(*i, useHitsInStat);
668 }
669 
670 static UInt_t maxTracks = 0;
671 static UInt_t maxMuchPts1 = 0;
672 static UInt_t maxMuchPts0 = 0;
673 static UInt_t maxStsPts7 = 0;
674 static UInt_t maxStsPts6 = 0;
675 
676 static inline void BuildStatistics(LxSimpleTrack* track) {
677  if (track->muchPoints[1][LXMIDDLE].size() > maxMuchPts1)
678  maxMuchPts1 = track->muchPoints[1][LXMIDDLE].size();
679 
680  if (track->muchPoints[0][LXMIDDLE].size() > maxMuchPts0)
681  maxMuchPts0 = track->muchPoints[0][LXMIDDLE].size();
682 
683  if (track->stsPoints[7].size() > maxStsPts7)
684  maxStsPts7 = track->stsPoints[7].size();
685 
686  if (track->stsPoints[6].size() > maxStsPts6)
687  maxStsPts6 = track->stsPoints[6].size();
688 
689  jPsiMuonsMomsHisto->Fill(track->p);
690 
691  for (list<LxSimplePoint>::iterator j = track->muchPoints[1][LXMIDDLE].begin();
692  j != track->muchPoints[1][LXMIDDLE].end();
693  ++j) {
694  LxSimplePoint muchPt1 = *j;
695 
696  for (list<LxSimplePoint>::iterator k =
697  track->muchPoints[0][LXMIDDLE].begin();
698  k != track->muchPoints[0][LXMIDDLE].end();
699  ++k) {
700  LxSimplePoint muchPt0 = *k;
701  Double_t diffZMuch = muchPt0.z - muchPt1.z;
702  Double_t txMuch = (muchPt0.x - muchPt1.x) / diffZMuch;
703  Double_t tyMuch = (muchPt0.y - muchPt1.y) / diffZMuch;
704 
705  if (-13 == track->pdgCode) {
706  Double_t txDiff = txMuch - muchPt0.x / muchPt0.z;
707  muPlusVertexTxDiff->Fill(txDiff);
708 
709  if (!track->stsPoints[0].empty() && !track->stsPoints[1].empty()) {
710  LxSimplePoint p0 = track->stsPoints[0].front();
711  LxSimplePoint p1 = track->stsPoints[1].front();
712  muPlusStsBeginTxDiff2D->Fill(track->p,
713  txMuch - (p1.x - p0.x) / (p1.z - p0.z));
714  muMinusStsBeginTxDiff2D->Fill(track->p,
715  (p1.x - p0.x) / (p1.z - p0.z) - txMuch);
716  deltaPhiPi->Fill(track->p * (txMuch - (p1.x - p0.x) / (p1.z - p0.z)));
717  }
718  } else if (13 == track->pdgCode) {
719  Double_t txDiff = txMuch - muchPt0.x / muchPt0.z;
720  muMinusVertexTxDiff->Fill(txDiff);
721 
722  if (!track->stsPoints[0].empty() && !track->stsPoints[1].empty()) {
723  LxSimplePoint p0 = track->stsPoints[0].front();
724  LxSimplePoint p1 = track->stsPoints[1].front();
725  muMinusStsBeginTxDiff2D->Fill(track->p,
726  txMuch - (p1.x - p0.x) / (p1.z - p0.z));
727  muPlusStsBeginTxDiff2D->Fill(track->p,
728  (p1.x - p0.x) / (p1.z - p0.z) - txMuch);
729  deltaPhiPi->Fill(track->p * ((p1.x - p0.x) / (p1.z - p0.z) - txMuch));
730  }
731  }
732 
733  for (list<LxSimplePoint>::iterator l = track->stsPoints[7].begin();
734  l != track->stsPoints[7].end();
735  ++l) {
736  LxSimplePoint stsPt7 = *l;
737  Double_t diffZ = stsPt7.z - muchPt0.z;
738  Double_t extX = muchPt0.x + txMuch * diffZ;
739  Double_t extY = muchPt0.y + tyMuch * diffZ;
740  Double_t dx = stsPt7.x - extX;
741  Double_t dy = stsPt7.y - extY;
742  muchStsBreakX->Fill(dx);
743  muchStsBreakY->Fill(dy);
744 
745  if (-13 == track->pdgCode) {
746  Double_t extXmu = muchPt0.x + (txMuch - 0.00671) * diffZ;
747  muPlusStsXDiff->Fill(stsPt7.x - extXmu);
748  } else if (13 == track->pdgCode) {
749  Double_t extXmu = muchPt0.x + (txMuch + 0.00691) * diffZ;
750  muMinusStsXDiff->Fill(stsPt7.x - extXmu);
751  }
752 
753  for (list<LxSimplePoint>::iterator m = track->stsPoints[6].begin();
754  m != track->stsPoints[6].end();
755  ++m) {
756  LxSimplePoint stsPt6 = *m;
757  Double_t diffZSts = stsPt6.z - stsPt7.z;
758  Double_t txSts = (stsPt6.x - stsPt7.x) / diffZSts;
759  Double_t tySts = (stsPt6.y - stsPt7.y) / diffZSts;
760  //muchStsBreakX->Fill(dx);
761  //muchStsBreakY->Fill(dy);
762  Double_t dtx = txSts - txMuch;
763  Double_t dty = tySts - tyMuch;
764  muchStsBreakTx->Fill(dtx);
765  muchStsBreakTy->Fill(dty);
766 
767  if (-13 == track->pdgCode)
768  muPlusStsTxDiff->Fill(dtx);
769  else if (13 == track->pdgCode)
770  muMinusStsTxDiff->Fill(dtx);
771 
772  Double_t chi2 = dx * dx / xRms2 + dy * dy / yRms2 + dtx * dtx / txRms2
773  + dty * dty / tyRms2;
774 
775  if (0 > track->motherId
776  && (13 == track->pdgCode || -13 == track->pdgCode)) // JPsi
777  signalChi2->Fill(chi2);
778  else
779  bgrChi2->Fill(chi2);
780 
781  diffZ = muchPt0.z - stsPt7.z;
782  extX = stsPt7.x + txSts * diffZ;
783  extY = stsPt7.y + tySts * diffZ;
784  dx = muchPt0.x - extX;
785  dy = muchPt0.y - extY;
786  stsMuchBreakX->Fill(dx);
787  stsMuchBreakY->Fill(dy);
788  }
789  }
790  }
791  }
792 }
793 
794 static inline void BuildNearestHitStat(LxSimpleTrack* track, bool cropHits) {
795  static Double_t maxDist = 0;
796  static Double_t mcX = 0;
797  static Double_t hitX = 0;
798  static Double_t mcY = 0;
799  static Double_t hitY = 0;
800 
801  static Double_t maxMinDist = 0;
802  static Double_t mcMinX = 0;
803  static Double_t hitMinX = 0;
804  static Double_t mcMinY = 0;
805  static Double_t hitMinY = 0;
806  static Int_t nMinHits = 0;
807 
808  for (Int_t i = 0; i < LXSTATIONS; ++i) {
809  if (track->muchMCPts[i][LXMIDDLE].empty()
810  || track->muchPoints[i][LXMIDDLE].empty())
811  continue;
812 
813  LxSimplePoint mcPoint =
814  track->muchMCPts[i][LXMIDDLE]
815  .front(); // We assume the MC points were averaged and there can be at most one point.
816  Double_t minDist = -1;
817  Double_t mcMinX0 = 0;
818  Double_t hitMinX0 = 0;
819  Double_t mcMinY0 = 0;
820  Double_t hitMinY0 = 0;
821  Double_t hitMinZ0 = 0;
822  Int_t nMinHits0 = 0;
823 
824  for (list<LxSimplePoint>::iterator j =
825  track->muchPoints[i][LXMIDDLE].begin();
826  j != track->muchPoints[i][LXMIDDLE].end();
827  ++j) {
828  LxSimplePoint hitPoint = *j;
829  Double_t dx = hitPoint.x - mcPoint.x;
830  Double_t dy = hitPoint.y - mcPoint.y;
831  Double_t dist = sqrt(dx * dx + dy * dy);
832 
833  if (track->motherId < 0
834  && (13 == track->pdgCode || -13 == track->pdgCode))
835  hitsDist[i]->Fill(dist);
836 
837  if (minDist > dist || minDist < 0) {
838  minDist = dist;
839 
840  mcMinX0 = mcPoint.x;
841  hitMinX0 = hitPoint.x;
842  mcMinY0 = mcPoint.y;
843  hitMinY0 = hitPoint.y;
844  hitMinZ0 = hitPoint.z;
845  nMinHits0 = track->muchPoints[i][LXMIDDLE].size();
846  }
847 
848  if (maxDist < dist) {
849  maxDist = dist;
850  mcX = mcPoint.x;
851  hitX = hitPoint.x;
852  mcY = mcPoint.y;
853  hitY = hitPoint.y;
854  }
855  } // for (list<LxSimplePoint>::iterator j = track->muchPoints[i].begin(); j != track->muchPoints[i].end(); ++j)
856 
857  if (minDist >= 0) {
858  if (track->motherId < 0
859  && (13 == track->pdgCode || -13 == track->pdgCode))
860  nearestHitDist[i]->Fill(minDist);
861 
862  if (cropHits) {
863  track->muchPoints[i][LXMIDDLE].clear();
864  LxSimplePoint point(hitMinX0, hitMinY0, hitMinZ0, 0, 0);
865  track->muchPoints[i][LXMIDDLE].push_back(point);
866  }
867 
868  if (maxMinDist < minDist) {
869  maxMinDist = minDist;
870  mcMinX = mcMinX0;
871  hitMinX = hitMinX0;
872  mcMinY = mcMinY0;
873  hitMinY = hitMinY0;
874  nMinHits = nMinHits0;
875  }
876  }
877  }
878 
879  cout << "BuildNearestHitStat: maxDist=" << maxDist << " MC x=" << mcX
880  << " Hit x=" << hitX << " MC y=" << mcY << " Hit y=" << hitY << endl;
881  cout << "BuildNearestHitStat: maxMinDist=" << maxMinDist << " MC x=" << mcMinX
882  << " Hit x=" << hitMinX << " MC y=" << mcMinY << " Hit y=" << hitMinY
883  << " n hits=" << nMinHits << endl;
884 }
885 
886 void LxTrackAnaTriplet::BuildStatistics() {
887  if (allTracks.size() > maxTracks) maxTracks = allTracks.size();
888 
889  for (vector<LxSimpleTrack*>::iterator i = allTracks.begin();
890  i != allTracks.end();
891  ++i) {
892  LxSimpleTrack* track = *i;
893 
894  if (track->muchPoints[0][LXMIDDLE].empty()
895  || track->muchPoints[1][LXMIDDLE].empty()
896  || track->muchPoints[5][LXMIDDLE].empty())
897  continue;
898 
899  if (buildNearestHitDist && useHitsInStat)
900  ::BuildNearestHitStat(
901  track,
902  cropHits); // Hits can be cropped (only the nearest to MC hit is left) here!
903 
904  if (track->motherId > -1 || (13 != track->pdgCode && -13 != track->pdgCode))
905  continue;
906 
907  ::BuildStatistics(track);
908  }
909 
910  cout << "Statistics is built maxTracks=" << maxTracks
911  << " maxMuchPts1=" << maxMuchPts1 << " maxMuchPts0=" << maxMuchPts0
912  << " maxStsPts7=" << maxStsPts7 << " maxStsPts6=" << maxStsPts6 << endl;
913 }
914 
915 void LxTrackAnaTriplet::Connect(bool useCuts) {
916  static Int_t jpsiTrackCount = 0;
917  static Int_t jpsiTrackCountCutted = 0;
918  static Int_t jpsiMatchedCount = 0;
919  static Int_t jpsiMatchedCountCutted = 0;
920 
921  static Int_t otherTrackCount = 0;
922  static Int_t otherTrackCountCutted = 0;
923  static Int_t otherMatchedCount = 0;
924  static Int_t otherMatchedCountCutted = 0;
925 
926  for (vector<LxSimpleTrack*>::iterator i = allTracks.begin();
927  i != allTracks.end();
928  ++i) {
929  LxSimpleTrack* muchTrack = *i;
930 
931  if (muchTrack->muchPoints[0][LXMIDDLE].empty()
932  || muchTrack->muchPoints[1][LXMIDDLE].empty()
933  || muchTrack->muchPoints[5][LXMIDDLE].empty())
934  continue;
935 
936  for (list<LxSimplePoint>::iterator j =
937  muchTrack->muchPoints[1][LXMIDDLE].begin();
938  j != muchTrack->muchPoints[1][LXMIDDLE].end();
939  ++j) {
940  LxSimplePoint muchPt1 = *j;
941 
942  for (list<LxSimplePoint>::iterator k =
943  muchTrack->muchPoints[0][LXMIDDLE].begin();
944  k != muchTrack->muchPoints[0][LXMIDDLE].end();
945  ++k) {
946  LxSimplePoint muchPt0 = *k;
947  Double_t diffZMuch = muchPt0.z - muchPt1.z;
948  Double_t txMuch = (muchPt0.x - muchPt1.x) / diffZMuch;
949  // Double_t txMuchVertex = muchPt0.x / muchPt0.z;
950  Double_t tyMuch = (muchPt0.y - muchPt1.y) / diffZMuch;
951  Connect(muchTrack, muchPt0, txMuch, tyMuch, useCuts);
952  } // for (list<LxSimplePoint>::iterator k = muchTrack->muchPoints[0].begin(); k != muchTrack->muchPoints[0].end(); ++k)
953  } // for (list<LxSimplePoint>::iterator j = muchTrack->muchPoints[1].begin(); j != muchTrack->muchPoints[1].end(); ++j)
954 
955  if (0 > muchTrack->motherId
956  && (13 == muchTrack->pdgCode || -13 == muchTrack->pdgCode)) // JPsi
957  {
958  if (useCuts) {
959  ++jpsiTrackCountCutted;
960 
961  if (muchTrack == muchTrack->linkedStsTrack) ++jpsiMatchedCountCutted;
962 
963  if (buildSigInvMass && 0 != muchTrack->linkedStsTrack) {
964  if (muchTrack->linkedStsTrack->charge > 0)
965  posTracks.push_back(muchTrack->linkedStsTrack);
966  else if (muchTrack->linkedStsTrack->charge < 0)
967  negTracks.push_back(muchTrack->linkedStsTrack);
968  }
969  } else {
970  ++jpsiTrackCount;
971 
972  if (muchTrack == muchTrack->linkedStsTrack) ++jpsiMatchedCount;
973  }
974  } else // Background track handled here.
975  {
976  if (useCuts) {
977  ++otherTrackCountCutted;
978 
979  if (muchTrack == muchTrack->linkedStsTrack
980  || 0 == muchTrack->linkedStsTrack)
981  ++otherMatchedCountCutted;
982 
983  if (buildBgrInvMass && !joinData && 0 != muchTrack->linkedStsTrack) {
984  superEventBrachTrack.px = muchTrack->linkedStsTrack->px;
985  superEventBrachTrack.py = muchTrack->linkedStsTrack->py;
986  superEventBrachTrack.pz = muchTrack->linkedStsTrack->pz;
987  superEventBrachTrack.e = muchTrack->linkedStsTrack->e;
988  superEventBrachTrack.charge = muchTrack->linkedStsTrack->charge;
989  superEventTracks->Fill();
990  }
991  } else {
992  ++otherTrackCount;
993 
994  if (muchTrack == muchTrack->linkedStsTrack
995  || 0 == muchTrack->linkedStsTrack)
996  ++otherMatchedCount;
997  }
998  }
999  } // for (vector<LxSimpleTrack*>::iterator i = allTracks.begin(); i != allTracks.end(); ++i)
1000 
1001  if (useCuts) {
1002  Double_t efficiency = jpsiMatchedCountCutted * 100;
1003  efficiency /= jpsiTrackCountCutted;
1004  cout << "J/psi (with cuts) connection efficiency = " << efficiency << "% ( "
1005  << jpsiMatchedCountCutted << " / " << jpsiTrackCountCutted << " )"
1006  << endl;
1007  efficiency = otherMatchedCountCutted * 100;
1008  efficiency /= otherTrackCountCutted;
1009  cout << "Others (with cuts) connection efficiency = " << efficiency
1010  << "% ( " << otherMatchedCountCutted << " / " << otherTrackCountCutted
1011  << " )" << endl;
1012  } else {
1013  Double_t efficiency = jpsiMatchedCount * 100;
1014  efficiency /= jpsiTrackCount;
1015  cout << "J/psi (without cuts) connection efficiency = " << efficiency
1016  << "% ( " << jpsiMatchedCount << " / " << jpsiTrackCount << " )"
1017  << endl;
1018  efficiency = otherMatchedCount * 100;
1019  efficiency /= otherTrackCount;
1020  cout << "Others (without cuts) connection efficiency = " << efficiency
1021  << "% ( " << otherMatchedCount << " / " << otherTrackCount << " )"
1022  << endl;
1023  }
1024 }
1025 
1026 void LxTrackAnaTriplet::Connect(LxSimpleTrack* muchTrack,
1027  LxSimplePoint muchPt0,
1028  Double_t txMuch,
1029  Double_t tyMuch,
1030  bool useCuts) {
1031  for (vector<LxSimpleTrack*>::iterator l = allTracks.begin();
1032  l != allTracks.end();
1033  ++l) {
1034  LxSimpleTrack* stsTrack = *l;
1035 
1036  if (stsTrack->p < 3.0 || stsTrack->pt < 1.0) continue;
1037 
1038  Int_t m0 = 0;
1039  Int_t n0 = -1;
1040 
1041  for (; m0 < LXSTSSTATIONS - 1; ++m0) {
1042  if (!stsTrack->stsPoints[m0].empty()) {
1043  n0 = m0 + 1;
1044 
1045  for (; n0 <= m0 + 2 && n0 < LXSTSSTATIONS; ++n0) {
1046  if (!stsTrack->stsPoints[n0].empty()) break;
1047  }
1048  }
1049 
1050  if (n0 <= m0 + 2) break;
1051  }
1052 
1053  Int_t m = LXSTSSTATIONS - 1;
1054  Int_t n = -1;
1055 
1056  for (; m > 0; --m) {
1057  if (!stsTrack->stsPoints[m].empty()) {
1058  n = m - 1;
1059 
1060  for (; n >= m - 2 && n >= 0; --n) {
1061  if (!stsTrack->stsPoints[n].empty()) break;
1062  }
1063  }
1064 
1065  if (n >= m - 2) break;
1066  }
1067 
1068  if (n >= 0 && m > n && n >= m - 2) {
1069  if (m0 >= LXSTSSTATIONS - 1 || n0 >= LXSTSSTATIONS || m0 >= n0) {
1070  m0 = n;
1071  n0 = m;
1072  }
1073 
1074  LxSimplePoint stsPtM0 = stsTrack->stsPoints[m0].front();
1075  LxSimplePoint stsPtN0 = stsTrack->stsPoints[n0].front();
1076  Double_t txSts0 = (stsPtN0.x - stsPtM0.x) / (stsPtN0.z - stsPtM0.z);
1077 
1078  for (list<LxSimplePoint>::iterator o = stsTrack->stsPoints[m].begin();
1079  o != stsTrack->stsPoints[m].end();
1080  ++o) {
1081  LxSimplePoint stsPtM = *o;
1082  Double_t deltaZ = stsPtM.z - muchPt0.z;
1083  Double_t extX = muchPt0.x + txMuch * deltaZ;
1084  Double_t extY = muchPt0.y + tyMuch * deltaZ;
1085  Double_t dx = stsPtM.x - extX;
1086  Double_t dy = stsPtM.y - extY;
1087 
1088  if (dx < 0) dx = -dx;
1089 
1090  if (dy < 0) dy = -dy;
1091 
1092  if (useCuts && (dx > xRms * cutCoeff || dy > yRms * cutCoeff)) continue;
1093 
1094  for (list<LxSimplePoint>::iterator p = stsTrack->stsPoints[n].begin();
1095  p != stsTrack->stsPoints[n].end();
1096  ++p) {
1097  LxSimplePoint stsPtN = *p;
1098  Double_t diffZSts = stsPtN.z - stsPtM.z;
1099  Double_t txSts = (stsPtN.x - stsPtM.x) / diffZSts;
1100  Double_t tySts = (stsPtN.y - stsPtM.y) / diffZSts;
1101  Double_t dtx = txSts - txMuch;
1102  Double_t dty = tySts - tyMuch;
1103 
1104  Double_t stsCharge =
1105  TDatabasePDG::Instance()->GetParticle(stsTrack->pdgCode)->Charge();
1106  Double_t muchCharge = txMuch - txSts0;
1107  bool chargesSignsTheSame = (stsCharge > 0 && muchCharge > 0)
1108  || (stsCharge < 0 && muchCharge < 0);
1109 
1110  if (dtx < 0) dtx = -dtx;
1111 
1112  if (dty < 0) dty = -dty;
1113 
1114  if (useCuts
1115  && ((useChargeSignInCuts
1116  && (!chargesSignsTheSame
1117  || !momFitTxBreak(stsTrack->p, muchCharge)))
1118  || dtx > txRms * cutCoeff || dty > tyRms * cutCoeff))
1119  continue;
1120 
1121  Double_t chi2 = dx * dx / xRms2 + dy * dy / yRms2 + dtx * dtx / txRms2
1122  + dty * dty / tyRms2;
1123  stsTrack->charge = stsCharge;
1124 
1125  if (0 == stsTrack->linkedMuchTrack.first
1126  || chi2 < stsTrack->linkedMuchTrack.second) {
1127  list<pair<LxSimpleTrack*, Double_t>>::iterator r =
1128  muchTrack->linkedStsTracks.begin();
1129 
1130  for (; r != muchTrack->linkedStsTracks.end() && r->second <= chi2;
1131  ++r)
1132  ;
1133 
1134  pair<LxSimpleTrack*, Double_t> trackDesc(stsTrack, chi2);
1135  muchTrack->linkedStsTracks.insert(r, trackDesc);
1136  }
1137  }
1138  }
1139  }
1140  } // for (vector<LxSimpleTrack*>::iterator l = allTracks.begin(); l != allTracks.end(); ++l)
1141 
1142  if (!muchTrack->linkedStsTracks.empty()) {
1143  pair<LxSimpleTrack*, Double_t> trackDesc =
1144  muchTrack->linkedStsTracks.front();
1145  muchTrack->linkedStsTracks.pop_front();
1146  LxSimpleTrack* anotherMuchTrack = trackDesc.first->linkedMuchTrack.first;
1147  trackDesc.first->linkedMuchTrack.first = muchTrack;
1148  trackDesc.first->linkedMuchTrack.second = trackDesc.second;
1149  muchTrack->linkedStsTrack = trackDesc.first;
1150 
1151  if (0 != anotherMuchTrack) anotherMuchTrack->RebindMuchTrack();
1152  }
1153 }
LxTrackAnaTriplet::cropHits
bool cropHits
Definition: riplet/LxTrackAna.h:134
yRms
static Double_t yRms
Definition: riplet/LxTrackAna.cxx:28
muchStsBreakTy
static TH1F * muchStsBreakTy
Definition: riplet/LxTrackAna.cxx:41
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
LxSimpleTrack::motherId
Int_t motherId
Definition: Simple/LxTrackAna.h:27
bgrInvMass
static TH1F * bgrInvMass
Definition: riplet/LxTrackAna.cxx:49
CbmTrack::GetParamLast
const FairTrackParam * GetParamLast() const
Definition: CbmTrack.h:62
CbmMuchPoint
Definition: CbmMuchPoint.h:21
CbmPixelHit::Position
void Position(TVector3 &pos) const
Copies hit position to pos.
Definition: CbmPixelHit.cxx:64
signalChi2
static TH1F * signalChi2
Definition: riplet/LxTrackAna.cxx:46
maxMuchPts1
static UInt_t maxMuchPts1
Definition: riplet/LxTrackAna.cxx:671
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
LxSimpleTrack::linkedMuchTrack
std::pair< LxSimpleTrack *, scaltype > linkedMuchTrack
Definition: Simple/LxTrackAna.h:59
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
MomVsTxRange::momHigh
Double_t momHigh
Definition: riplet/LxTrackAna.cxx:73
stsMuchBreakX
static TH1F * stsMuchBreakX
Definition: riplet/LxTrackAna.cxx:43
BuildInvMass2
static void BuildInvMass2(list< CbmStsTrack * > &stsTracks, TH1 *)
Definition: riplet/LxTrackAna.cxx:335
LxSimpleTrack::py
scaltype py
Definition: Simple/LxTrackAna.h:31
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
muPlusStsTxDiff
static TH1F * muPlusStsTxDiff
Definition: riplet/LxTrackAna.cxx:57
muchStsBreakX
static TH1F * muchStsBreakX
Definition: riplet/LxTrackAna.cxx:38
CbmMuchCluster
Data container for MUCH clusters.
Definition: CbmMuchCluster.h:20
LxSimpleTrack::linkedStsTracks
std::list< std::pair< LxSimpleTrack *, scaltype > > linkedStsTracks
Definition: Simple/LxTrackAna.h:61
saveHistos
bool saveHistos
Definition: riplet/LxTrackAnaSegments.cxx:11
LxTrackAnaTriplet::Exec
void Exec(Option_t *opt)
Definition: riplet/LxTrackAna.cxx:472
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
LxTrackAnaTriplet::buildBgrInvMass
bool buildBgrInvMass
Definition: riplet/LxTrackAna.h:130
LxTrackAnaTriplet::segmentsAnalyzer
LxTrackAnaSegments segmentsAnalyzer
Definition: riplet/LxTrackAna.h:137
CbmMuchPoint::GetDetectorId
Int_t GetDetectorId() const
Definition: CbmMuchPoint.h:69
LxSimpleTrack::pz
scaltype pz
Definition: Simple/LxTrackAna.h:32
LxTrackAnaTriplet::Connect
void Connect(bool useCuts)
Definition: riplet/LxTrackAna.cxx:915
LxSimpleTrack::RebindMuchTrack
void RebindMuchTrack()
Definition: Simple/LxTrackAna.cxx:95
LxTrackAnaTriplet::listMuchPts
TClonesArray * listMuchPts
Definition: riplet/LxTrackAna.h:116
BuildNearestHitStat
static void BuildNearestHitStat(LxSimpleTrack *track, bool cropHits)
Definition: riplet/LxTrackAna.cxx:794
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
LXLAYERS
#define LXLAYERS
Definition: Simple/LxSettings.h:8
LxTrackAnaTriplet::listMuchPixelHits
TClonesArray * listMuchPixelHits
Definition: riplet/LxTrackAna.h:117
CbmMCTrack::GetPx
Double_t GetPx() const
Definition: CbmMCTrack.h:72
LxSimplePoint::y
scaltype y
Definition: Simple/LxTrackAna.h:16
LxSimpleTrack::pt
scaltype pt
Definition: Simple/LxTrackAna.h:29
CbmMCTrack::GetPy
Double_t GetPy() const
Definition: CbmMCTrack.h:73
LxTrackAnaTriplet::negTracks
std::list< LxSimpleTrack * > negTracks
Definition: riplet/LxTrackAna.h:122
MomVsTxRange
Definition: Simple/LxTrackAna.cxx:79
CbmStsPoint
Definition: CbmStsPoint.h:27
xRms
static Double_t xRms
Definition: riplet/LxTrackAna.cxx:25
LXMIDDLE
#define LXMIDDLE
Definition: Simple/LxSettings.h:10
LxSimpleTrack
Definition: Simple/LxTrackAna.h:25
LxTrackAnaTriplet
Definition: riplet/LxTrackAna.h:68
LxTrackAnaTriplet::joinData
bool joinData
Definition: riplet/LxTrackAna.h:132
yRms2
static Double_t yRms2
Definition: riplet/LxTrackAna.cxx:29
LxTrackAnaTriplet::FinishTask
void FinishTask()
Definition: riplet/LxTrackAna.cxx:381
CbmMuchDigiMatch
Definition: CbmMuchDigiMatch.h:17
CbmCluster::GetNofDigis
Int_t GetNofDigis() const
Number of digis in cluster.
Definition: CbmCluster.h:69
nearestHitDist
static TH1F * nearestHitDist[LXSTATIONS]
Definition: riplet/LxTrackAna.cxx:54
CbmMuchCluster.h
Data container for MUCH clusters.
LxTrackAna.h
particleType
static TString particleType("jpsi")
tyRms2
static Double_t tyRms2
Definition: riplet/LxTrackAna.cxx:35
CbmStsPoint::PositionOut
void PositionOut(TVector3 &pos)
Definition: CbmStsPoint.h:96
maxMuchPts0
static UInt_t maxMuchPts0
Definition: riplet/LxTrackAna.cxx:672
CbmStsTrack.h
Data class for STS tracks.
tyRms
static Double_t tyRms
Definition: riplet/LxTrackAna.cxx:34
CbmMuchGeoScheme::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchGeoScheme.h:71
LxTrackAnaTriplet::Init
InitStatus Init()
Definition: riplet/LxTrackAna.cxx:147
CbmMuchGeoScheme::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchGeoScheme.h:68
maxTracks
static UInt_t maxTracks
Definition: riplet/LxTrackAna.cxx:670
CbmMuchPoint.h
BuildStatistics
static void BuildStatistics(LxSimpleTrack *track)
Definition: riplet/LxTrackAna.cxx:676
LxSimplePoint::z
scaltype z
Definition: Simple/LxTrackAna.h:17
muMinusStsXDiff
static TH1F * muMinusStsXDiff
Definition: riplet/LxTrackAna.cxx:60
tracks
TClonesArray * tracks
Definition: Analyze_matching.h:17
muchStsBreakTx
static TH1F * muchStsBreakTx
Definition: riplet/LxTrackAna.cxx:40
LxTrackAnaSegments::BuildStatistics
void BuildStatistics()
Definition: Simple/LxTrackAnaSegments.cxx:353
LXSTSSTATIONS
#define LXSTSSTATIONS
Definition: Simple/LxSettings.h:12
CbmStsAddress::GetElementId
UInt_t GetElementId(Int_t address, Int_t level)
Get the index of an element.
Definition: CbmStsAddress.cxx:180
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
LxSimpleTrack::stsPoints
std::list< LxSimplePoint > stsPoints[LXSTSSTATIONS]
Definition: Simple/LxTrackAna.h:54
CbmStsPoint::MomentumOut
void MomentumOut(TVector3 &mom)
Definition: CbmStsPoint.h:97
LxTrackAnaTriplet::~LxTrackAnaTriplet
~LxTrackAnaTriplet()
Definition: riplet/LxTrackAna.cxx:134
LxTrackAnaTriplet::buildSegmentsStat
bool buildSegmentsStat
Definition: riplet/LxTrackAna.h:135
LxSimplePoint::x
scaltype x
Definition: Simple/LxTrackAna.h:15
txRms2
static Double_t txRms2
Definition: riplet/LxTrackAna.cxx:32
SaveHisto
static void SaveHisto(TH1 *histo, const char *particleType, const char *name)
Definition: riplet/LxTrackAna.cxx:294
CbmMuchPoint::PositionOut
void PositionOut(TVector3 &pos) const
Definition: CbmMuchPoint.h:80
CbmKFTrack.h
LxSimpleTrack::muchMCPts
std::list< LxSimplePoint > muchMCPts[LXSTATIONS][LXLAYERS]
Definition: Simple/LxTrackAna.h:58
LxTrackAnaTriplet::averagePoints
bool averagePoints
Definition: riplet/LxTrackAna.h:126
LxTrackAnaTriplet::useChargeSignInCuts
bool useChargeSignInCuts
Definition: riplet/LxTrackAna.h:128
LxSimpleTrack::charge
scaltype charge
Definition: Simple/LxTrackAna.h:34
hitsDist
static TH1F * hitsDist[LXSTATIONS]
Definition: riplet/LxTrackAna.cxx:55
CbmMuchPoint::MomentumOut
void MomentumOut(TVector3 &mom) const
Definition: CbmMuchPoint.h:81
LxTrackAnaTriplet::buildConnectStat
bool buildConnectStat
Definition: riplet/LxTrackAna.h:129
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
momFitTxBreak
static bool momFitTxBreak(Double_t mom, Double_t txBreak)
Definition: riplet/LxTrackAna.cxx:78
LxTrackAnaTriplet::LxTrackAnaTriplet
LxTrackAnaTriplet()
Definition: riplet/LxTrackAna.cxx:108
CbmKFParticle.h
muPlusVertexTxDiff
static TH1F * muPlusVertexTxDiff
Definition: riplet/LxTrackAna.cxx:61
CbmMCTrack::GetPt
Double_t GetPt() const
Definition: CbmMCTrack.h:99
muMinusStsTxDiff
static TH1F * muMinusStsTxDiff
Definition: riplet/LxTrackAna.cxx:58
MomVsTxRange::momLow
Double_t momLow
Definition: riplet/LxTrackAna.cxx:72
LxTrackAnaTriplet::superEventTracks
TTree * superEventTracks
Definition: riplet/LxTrackAna.h:123
ClassImp
ClassImp(LxTrackAnaTriplet) using namespace std
muMinusVertexTxDiff
static TH1F * muMinusVertexTxDiff
Definition: riplet/LxTrackAna.cxx:62
LxTrackAnaSegments::Init
void Init()
Definition: Simple/LxTrackAnaSegments.cxx:62
LxTrackAnaTriplet::Clean
void Clean()
Definition: riplet/LxTrackAna.cxx:136
bgrChi2
static TH1F * bgrChi2
Definition: riplet/LxTrackAna.cxx:47
CbmStsPoint.h
points
TClonesArray * points
Definition: Analyze_matching.h:18
CbmMCTrack.h
LxTrackAnaTriplet::AveragePoints
void AveragePoints()
Definition: riplet/LxTrackAna.cxx:663
LxSimpleTrack::parent
LxSimpleTrack * parent
Definition: riplet/LxTrackAna.h:64
muPlusStsBeginTxDiff2D
static TH2F * muPlusStsBeginTxDiff2D
Definition: riplet/LxTrackAna.cxx:64
LxTrackAnaTriplet::useHitsInStat
bool useHitsInStat
Definition: riplet/LxTrackAna.h:125
xRms2
static Double_t xRms2
Definition: riplet/LxTrackAna.cxx:26
LxSimpleTrack::e
scaltype e
Definition: Simple/LxTrackAna.h:33
negativeTracks
static list< LxSimpleTrack * > negativeTracks
Definition: riplet/LxTrackAna.cxx:51
CbmMCTrack
Definition: CbmMCTrack.h:34
LxSimplePoint::ty
scaltype ty
Definition: Simple/LxTrackAna.h:19
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LxSimpleTrack::linkedStsTrack
LxSimpleTrack * linkedStsTrack
Definition: Simple/LxTrackAna.h:62
txRms
static Double_t txRms
Definition: riplet/LxTrackAna.cxx:31
LxTrackAnaTriplet::particleType
TString particleType
Definition: riplet/LxTrackAna.h:136
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
LxTrackAnaTriplet::listStsPts
TClonesArray * listStsPts
Definition: riplet/LxTrackAna.h:115
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LxTrackAnaTriplet::listMuchPixelDigiMatches
TClonesArray * listMuchPixelDigiMatches
Definition: riplet/LxTrackAna.h:119
LxSimpleTrack::muchPoints
std::list< LxSimplePoint > muchPoints[LXSTATIONS][LXLAYERS]
Definition: Simple/LxTrackAna.h:55
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmStsAddress.h
BuildInvMass
static void BuildInvMass(list< LxSimpleTrack * > &pTracks, list< LxSimpleTrack * > &nTracks, TH1 *histo)
Definition: riplet/LxTrackAna.cxx:314
LxTrackAnaTriplet::listMuchClusters
TClonesArray * listMuchClusters
Definition: riplet/LxTrackAna.h:118
cutCoeff
static Double_t cutCoeff
Definition: riplet/LxTrackAna.cxx:36
CbmMuchDigiMatch.h
jPsiMuonsMomsHisto
static TH1F * jPsiMuonsMomsHisto
Definition: riplet/LxTrackAna.cxx:69
LxSimplePoint::tx
scaltype tx
Definition: Simple/LxTrackAna.h:18
LxTrackAnaTriplet::superEventBrachTrack
LxSimpleTrack superEventBrachTrack
Definition: riplet/LxTrackAna.h:124
CbmMuchPixelHit
Definition: CbmMuchPixelHit.h:17
CbmMuchGeoScheme.h
muchStsBreakY
static TH1F * muchStsBreakY
Definition: riplet/LxTrackAna.cxx:39
CbmMCTrack::GetEnergy
Double_t GetEnergy() const
Definition: CbmMCTrack.h:165
LxSimpleTrack::px
scaltype px
Definition: Simple/LxTrackAna.h:30
LxTrackAnaTriplet::posTracks
std::list< LxSimpleTrack * > posTracks
Definition: riplet/LxTrackAna.h:121
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmKFTrack
Definition: CbmKFTrack.h:21
LxSimpleTrack::p
scaltype p
Definition: Simple/LxTrackAna.h:28
positiveTracks
static list< LxSimpleTrack * > positiveTracks
Definition: riplet/LxTrackAna.cxx:50
LxTrackAnaTriplet::buildNearestHitDist
bool buildNearestHitDist
Definition: riplet/LxTrackAna.h:133
LxTrackAnaSegments::Finish
void Finish()
Definition: Simple/LxTrackAnaSegments.cxx:307
LXSTATIONS
#define LXSTATIONS
Definition: Simple/LxSettings.h:9
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
MomVsTxRange::txLow
Double_t txLow
Definition: riplet/LxTrackAna.cxx:74
muMinusStsBeginTxDiff2D
static TH2F * muMinusStsBeginTxDiff2D
Definition: riplet/LxTrackAna.cxx:65
CbmPixelHit::PositionError
void PositionError(TVector3 &dpos) const
Copies hit position error to pos.
Definition: CbmPixelHit.cxx:68
LxTrackAnaTriplet::listMCTracks
TClonesArray * listMCTracks
Definition: riplet/LxTrackAna.h:114
LxSimpleTrack::pdgCode
Int_t pdgCode
Definition: Simple/LxTrackAna.h:26
CbmMCTrack::GetP
Double_t GetP() const
Definition: CbmMCTrack.h:100
MomVsTxRange::txHigh
Double_t txHigh
Definition: riplet/LxTrackAna.cxx:75
sigInvMass
static TH1F * sigInvMass
Definition: riplet/LxTrackAna.cxx:52
muPlusStsXDiff
static TH1F * muPlusStsXDiff
Definition: riplet/LxTrackAna.cxx:59
stsMuchBreakY
static TH1F * stsMuchBreakY
Definition: riplet/LxTrackAna.cxx:44
LxTrackAnaTriplet::buildSigInvMass
bool buildSigInvMass
Definition: riplet/LxTrackAna.h:131
maxStsPts7
static UInt_t maxStsPts7
Definition: riplet/LxTrackAna.cxx:673
CbmMCTrack::GetPz
Double_t GetPz() const
Definition: CbmMCTrack.h:74
deltaPhiPi
static TH1F * deltaPhiPi
Definition: riplet/LxTrackAna.cxx:67
AveragePoints
static void AveragePoints(list< LxSimplePoint > &points)
Definition: riplet/LxTrackAna.cxx:616
LxTrackAnaTriplet::allTracks
std::vector< LxSimpleTrack * > allTracks
Definition: riplet/LxTrackAna.h:120
LxSimplePoint
Definition: Simple/LxTrackAna.h:14
maxStsPts6
static UInt_t maxStsPts6
Definition: riplet/LxTrackAna.cxx:674
LxTrackAnaTriplet::BuildStatistics
void BuildStatistics()
Definition: riplet/LxTrackAna.cxx:886