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