CbmRoot
riplet/LxTrackAnaSegments.cxx
Go to the documentation of this file.
1 #include "LxTrackAnaSegments.h"
2 #include "LxTrackAna.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include <cmath>
6 #include <dirent.h>
7 #include <iostream>
8 #include <sys/stat.h>
9 #include <sys/types.h>
10 
11 bool saveHistos = true;
12 
13 using namespace std;
14 
15 static TH1F* muchInStationXDispLeft[LXSTATIONS];
16 static TH1F* muchInStationXDispRight[LXSTATIONS];
17 static TH1F* muchInStationYDispLeft[LXSTATIONS];
18 static TH1F* muchInStationYDispRight[LXSTATIONS];
19 
20 static TH1F* muchInStationXDispRL[LXSTATIONS];
21 static TH1F* muchInStationYDispRL[LXSTATIONS];
22 
23 static TH1F* muchInStationTxBreak[LXSTATIONS];
24 static TH1F* muchInStationTyBreak[LXSTATIONS];
25 
26 static TH1F* muchOutStationTxBreakLeft[LXSTATIONS - 1];
27 static TH1F* muchOutStationTxBreakRight[LXSTATIONS - 1];
28 static TH1F* muchOutStationTyBreakLeft[LXSTATIONS - 1];
29 static TH1F* muchOutStationTyBreakRight[LXSTATIONS - 1];
30 
31 static TH1F* muchTripletTxBreak[LXSTATIONS - 1];
32 static TH1F* muchTripletTyBreak[LXSTATIONS - 1];
33 
34 static TH1F* muchOutStationXDispByTriplet[LXSTATIONS - 1];
35 static TH1F* muchOutStationYDispByTriplet[LXSTATIONS - 1];
36 
37 static TH1F* muchOutStationXDispByVertex[LXSTATIONS - 1];
38 static TH1F* muchOutStationYDispByVertex[LXSTATIONS - 1];
39 
40 static TH1F* muchLongSegmentTxHisto[LXSTATIONS - 1];
41 static TH1F* muchLongSegmentTyHisto[LXSTATIONS - 1];
42 
43 static TH1F* muchSegmentTxBreakHisto[LXSTATIONS - 2];
44 static TH1F* muchSegmentTyBreakHisto[LXSTATIONS - 2];
45 
46 static TH1F* muchStationTxDispHisto[LXSTATIONS - 1];
47 static TH1F* muchStationTyDispHisto[LXSTATIONS - 1];
48 
49 static TH2F* muchXTxCovHisto[LXSTATIONS - 1];
50 static TH2F* muchYTyCovHisto[LXSTATIONS - 1];
51 
52 static TH1F* muchClusterXDispHisto[LXSTATIONS - 1];
53 static TH1F* muchClusterYDispHisto[LXSTATIONS - 1];
54 static TH1F* muchClusterTxDispHisto[LXSTATIONS - 1];
55 static TH1F* muchClusterTyDispHisto[LXSTATIONS - 1];
56 
57 static bool
58 GetHistoRMS(const char* histoNameBase, Int_t histoNumber, Double_t& retVal) {
59  char name[256];
60  char dir_name[256];
61  sprintf(dir_name, "configuration.%s", "omega");
62  bool result = false;
63  sprintf(name, "%s/%s_%d.root", dir_name, histoNameBase, histoNumber);
64  TFile* curFile = TFile::CurrentFile();
65  TFile* f = new TFile(name);
66 
67  if (!f->IsZombie()) {
68  sprintf(name, "%s_%d", histoNameBase, histoNumber);
69  TH1F* h = static_cast<TH1F*>(f->Get(name));
70  retVal = h->GetRMS();
71  result = true;
72  }
73 
74  delete f;
75  TFile::CurrentFile() = curFile;
76  return result;
77 }
78 
79 LxTrackAnaSegments::LxTrackAnaSegments(LxTrackAnaTriplet& o)
80  : owner(o), stationsInAlgo(LXSTATIONS), useBgr(false) {}
81 
82 static TString particleType("jpsi");
83 
84 void LxTrackAnaSegments::SetParticleType(TString v) {
85  particleType = v;
86 
87  if (v == "omega") stationsInAlgo = 5;
88 }
89 
90 void LxTrackAnaSegments::Init() {
91  char name[64];
92  char title[256];
93 
94  for (Int_t i = 0; i < LXSTATIONS; ++i) {
95  sprintf(name, "muchInStationXDispLeft_%d", i);
96  sprintf(
97  title, "X dispersion from central to left layer inside station: %d", i);
98  muchInStationXDispLeft[i] = new TH1F(name, title, 100, -3.0, 3.0);
99  muchInStationXDispLeft[i]->StatOverflows();
100 
101  sprintf(name, "muchInStationXDispRight_%d", i);
102  sprintf(
103  title, "X dispersion from central to right layer inside station: %d", i);
104  muchInStationXDispRight[i] = new TH1F(name, title, 100, -3.0, 3.0);
105  muchInStationXDispRight[i]->StatOverflows();
106 
107  sprintf(name, "muchInStationYDispLeft_%d", i);
108  sprintf(
109  title, "Y dispersion from central to left layer inside station: %d", i);
110  muchInStationYDispLeft[i] = new TH1F(name, title, 100, -3.0, 3.0);
111  muchInStationYDispLeft[i]->StatOverflows();
112 
113  sprintf(name, "muchInStationYDispRight_%d", i);
114  sprintf(
115  title, "Y dispersion from central to right layer inside station: %d", i);
116  muchInStationYDispRight[i] = new TH1F(name, title, 100, -3.0, 3.0);
117  muchInStationYDispRight[i]->StatOverflows();
118 
119  sprintf(name, "muchInStationXDispRL_%d", i);
120  sprintf(
121  title, "X dispersion on left layer predicted by right station: %d", i);
122  muchInStationXDispRL[i] = new TH1F(name, title, 100, -0.1, 0.1);
123  muchInStationXDispRL[i]->StatOverflows();
124 
125  sprintf(name, "muchInStationYDispRL_%d", i);
126  sprintf(
127  title, "Y dispersion on left layer predicted by right station: %d", i);
128  muchInStationYDispRL[i] = new TH1F(name, title, 100, -0.1, 0.1);
129  muchInStationYDispRL[i]->StatOverflows();
130 
131  sprintf(name, "muchInStationTxBreak_%d", i);
132  sprintf(title, "Tx break inside station: %d", i);
133  muchInStationTxBreak[i] = new TH1F(name, title, 100, -0.02, 0.02);
134  muchInStationTxBreak[i]->StatOverflows();
135 
136  sprintf(name, "muchInStationTyBreak_%d", i);
137  sprintf(title, "Ty break inside station: %d", i);
138  muchInStationTyBreak[i] = new TH1F(name, title, 100, -0.02, 0.02);
139  muchInStationTyBreak[i]->StatOverflows();
140 
141  if (i > 0) {
142  sprintf(name, "muchLongSegmentTxHisto_%d", i);
143  sprintf(
144  title,
145  "Tx tangents distribution for segments between stations: %d and %d",
146  i - 1,
147  i);
148  muchLongSegmentTxHisto[i - 1] = new TH1F(name, title, 100, -.15, .15);
149  muchLongSegmentTxHisto[i - 1]->StatOverflows();
150 
151  sprintf(name, "muchLongSegmentTyHisto_%d", i);
152  sprintf(
153  title,
154  "Ty tangents distribution for segments between stations: %d and %d",
155  i - 1,
156  i);
157  muchLongSegmentTyHisto[i - 1] = new TH1F(name, title, 100, -.15, .15);
158  muchLongSegmentTyHisto[i - 1]->StatOverflows();
159 
160  sprintf(name, "muchClusterXDispHisto_%d", i);
161  sprintf(title,
162  "X coordinate dispersion for cluster segments between stations: "
163  "%d and %d",
164  i - 1,
165  i);
166  muchClusterXDispHisto[i - 1] = new TH1F(name, title, 100, .0, 3.0);
167  muchClusterXDispHisto[i - 1]->StatOverflows();
168 
169  sprintf(name, "muchClusterYDispHisto_%d", i);
170  sprintf(title,
171  "Y coordinate dispersion for cluster segments between stations: "
172  "%d and %d",
173  i - 1,
174  i);
175  muchClusterYDispHisto[i - 1] = new TH1F(name, title, 100, .0, 3.0);
176  muchClusterYDispHisto[i - 1]->StatOverflows();
177 
178  sprintf(name, "muchClusterTxDispHisto_%d", i);
179  sprintf(title,
180  "Tx tangent dispersion for cluster segments between stations: %d "
181  "and %d",
182  i - 1,
183  i);
184  muchClusterTxDispHisto[i - 1] = new TH1F(name, title, 100, .0, .05);
185  muchClusterTxDispHisto[i - 1]->StatOverflows();
186 
187  sprintf(name, "muchClusterTyDispHisto_%d", i);
188  sprintf(title,
189  "Ty tangent dispersion for cluster segments between stations: %d "
190  "and %d",
191  i - 1,
192  i);
193  muchClusterTyDispHisto[i - 1] = new TH1F(name, title, 100, .0, .05);
194  muchClusterTyDispHisto[i - 1]->StatOverflows();
195 
196  sprintf(name, "muchOutStationTxBreakLeft_%d", i);
197  sprintf(title,
198  "Tx break between right segment of station and left tip of the "
199  "interstation segment: %d",
200  i);
201  muchOutStationTxBreakLeft[i - 1] =
202  new TH1F(name, title, 100, -0.15, 0.15);
203  muchOutStationTxBreakLeft[i - 1]->StatOverflows();
204 
205  sprintf(name, "muchOutStationTxBreakRight_%d", i);
206  sprintf(title,
207  "Tx break between left segment of station and right tip of the "
208  "interstation segment: %d",
209  i);
210  muchOutStationTxBreakRight[i - 1] =
211  new TH1F(name, title, 100, -0.15, 0.15);
212  muchOutStationTxBreakRight[i - 1]->StatOverflows();
213 
214  sprintf(name, "muchOutStationTyBreakLeft_%d", i);
215  sprintf(title,
216  "Ty break between right segment of station and left tip of the "
217  "interstation segment: %d",
218  i);
219  muchOutStationTyBreakLeft[i - 1] =
220  new TH1F(name, title, 100, -0.15, 0.15);
221  muchOutStationTyBreakLeft[i - 1]->StatOverflows();
222 
223  sprintf(name, "muchOutStationTyBreakRight_%d", i);
224  sprintf(title,
225  "Ty break between left segment of station and right tip of the "
226  "interstation segment: %d",
227  i);
228  muchOutStationTyBreakRight[i - 1] =
229  new TH1F(name, title, 100, -0.15, 0.15);
230  muchOutStationTyBreakRight[i - 1]->StatOverflows();
231 
232  sprintf(name, "muchOutStationXDispByTriplet_%d", i);
233  sprintf(title,
234  "X dispersion of prediction by triplet angle for station: %d",
235  i);
236  muchOutStationXDispByTriplet[i - 1] =
237  new TH1F(name, title, 100, -10.0, 10.0);
238  muchOutStationXDispByTriplet[i - 1]->StatOverflows();
239 
240  sprintf(name, "muchOutStationYDispByTriplet_%d", i);
241  sprintf(title,
242  "Y dispersion of prediction by triplet angle for station: %d",
243  i);
244  muchOutStationYDispByTriplet[i - 1] =
245  new TH1F(name, title, 100, -10.0, 10.0);
246  muchOutStationYDispByTriplet[i - 1]->StatOverflows();
247 
248  sprintf(name, "muchOutStationXDispByVertex_%d", i);
249  sprintf(
250  title,
251  "X dispersion of prediction by an angle to vertex for station: %d",
252  i);
253  muchOutStationXDispByVertex[i - 1] =
254  new TH1F(name, title, 100, -10.0, 10.0);
255  muchOutStationXDispByVertex[i - 1]->StatOverflows();
256 
257  sprintf(name, "muchOutStationYDispByVertex_%d", i);
258  sprintf(
259  title,
260  "Y dispersion of prediction by an angle to vertex for station: %d",
261  i);
262  muchOutStationYDispByVertex[i - 1] =
263  new TH1F(name, title, 100, -10.0, 10.0);
264  muchOutStationYDispByVertex[i - 1]->StatOverflows();
265 
266  sprintf(name, "muchTripletTxBreak_%d", i);
267  sprintf(
268  title, "Tx break between triplets on stations %d and %d", i - 1, i);
269  muchTripletTxBreak[i - 1] = new TH1F(name, title, 200, -0.2, 0.2);
270  muchTripletTxBreak[i - 1]->StatOverflows();
271 
272  sprintf(name, "muchTripletTyBreak_%d", i);
273  sprintf(
274  title, "Ty break between triplets on stations %d and %d", i - 1, i);
275  muchTripletTyBreak[i - 1] = new TH1F(name, title, 200, -0.2, 0.2);
276  muchTripletTyBreak[i - 1]->StatOverflows();
277 
278  if (i < LXSTATIONS - 1) {
279  sprintf(name, "muchSegmentTxBreakHisto_%d", i);
280  sprintf(title,
281  "Tx tangents breaks distribution for adjacent segments on "
282  "station: %d",
283  i);
284  muchSegmentTxBreakHisto[i - 1] = new TH1F(name, title, 100, -.15, .15);
285  muchSegmentTxBreakHisto[i - 1]->StatOverflows();
286 
287  sprintf(name, "muchSegmentTyBreakHisto_%d", i);
288  sprintf(title,
289  "Ty tangents breaks distribution for adjacent segments on "
290  "station: %d",
291  i);
292  muchSegmentTyBreakHisto[i - 1] = new TH1F(name, title, 100, -.15, .15);
293  muchSegmentTyBreakHisto[i - 1]->StatOverflows();
294  }
295 
296  sprintf(name, "muchStationTxDispHisto_%d", i);
297  sprintf(title,
298  "Tx tangents dispersion for segments between stations: %d and %d",
299  i - 1,
300  i);
301  muchStationTxDispHisto[i - 1] = new TH1F(name, title, 100, -.05, .05);
302  muchStationTxDispHisto[i - 1]->StatOverflows();
303 
304  sprintf(name, "muchStationTyDispHisto_%d", i);
305  sprintf(title,
306  "Ty tangents dispersion for segments between stations: %d and %d",
307  i - 1,
308  i);
309  muchStationTyDispHisto[i - 1] = new TH1F(name, title, 100, -.05, .05);
310  muchStationTyDispHisto[i - 1]->StatOverflows();
311  }
312 
313  if (i < LXSTATIONS - 1) {
314  sprintf(name, "muchXTxCovHisto_%d", i);
315  sprintf(title, "muchXTxCovHisto on %d", i);
316  muchXTxCovHisto[i] =
317  new TH2F(name, title, 100, -5.0, 5.0, 100, -0.15, 0.15);
318  muchXTxCovHisto[i]->StatOverflows();
319 
320  sprintf(name, "muchYTyCovHisto_%d", i);
321  sprintf(title, "muchYTyCovHisto on %d", i);
322  muchYTyCovHisto[i] =
323  new TH2F(name, title, 100, -5.0, 5.0, 100, -0.15, 0.15);
324  muchYTyCovHisto[i]->StatOverflows();
325  }
326  }
327 }
328 
329 static void SaveHisto(TH1* histo) {
330  if (!saveHistos) return;
331 
332  char dir_name[256];
333  sprintf(dir_name, "configuration.%s", particleType.Data());
334  DIR* dir = opendir(dir_name);
335 
336  if (dir)
337  closedir(dir);
338  else
339  mkdir(dir_name, 0700);
340 
341  char name[256];
342  sprintf(name, "%s/%s.root", dir_name, histo->GetName());
343  TFile fh(name, "RECREATE");
344  histo->Write();
345  fh.Close();
346  delete histo;
347 }
348 
349 void LxTrackAnaSegments::Finish() {
350  for (Int_t i = 0; i < LXSTATIONS; ++i) {
351  SaveHisto(muchInStationXDispLeft[i]);
352  SaveHisto(muchInStationXDispRight[i]);
353  SaveHisto(muchInStationYDispLeft[i]);
354  SaveHisto(muchInStationYDispRight[i]);
355  SaveHisto(muchInStationXDispRL[i]);
356  SaveHisto(muchInStationYDispRL[i]);
357  SaveHisto(muchInStationTxBreak[i]);
358  SaveHisto(muchInStationTyBreak[i]);
359 
360  if (i > 0) {
361  SaveHisto(muchLongSegmentTxHisto[i - 1]);
362  SaveHisto(muchLongSegmentTyHisto[i - 1]);
363 
364  SaveHisto(muchClusterXDispHisto[i - 1]);
365  SaveHisto(muchClusterYDispHisto[i - 1]);
366  SaveHisto(muchClusterTxDispHisto[i - 1]);
367  SaveHisto(muchClusterTyDispHisto[i - 1]);
368 
369  SaveHisto(muchOutStationTxBreakLeft[i - 1]);
370  SaveHisto(muchOutStationTxBreakRight[i - 1]);
371  SaveHisto(muchOutStationTyBreakLeft[i - 1]);
372  SaveHisto(muchOutStationTyBreakRight[i - 1]);
373 
374  SaveHisto(muchOutStationXDispByTriplet[i - 1]);
375  SaveHisto(muchOutStationYDispByTriplet[i - 1]);
376  SaveHisto(muchOutStationXDispByVertex[i - 1]);
377  SaveHisto(muchOutStationYDispByVertex[i - 1]);
378 
379  SaveHisto(muchTripletTxBreak[i - 1]);
380  SaveHisto(muchTripletTyBreak[i - 1]);
381 
382  if (i < LXSTATIONS - 1) {
383  SaveHisto(muchSegmentTxBreakHisto[i - 1]);
384  SaveHisto(muchSegmentTyBreakHisto[i - 1]);
385  }
386 
387  SaveHisto(muchStationTxDispHisto[i - 1]);
388  SaveHisto(muchStationTyDispHisto[i - 1]);
389  }
390 
391  if (i < LXSTATIONS - 1) {
392  SaveHisto(muchXTxCovHisto[i]);
393  SaveHisto(muchYTyCovHisto[i]);
394  }
395  }
396 }
397 
398 static bool GetPoints(LxSimpleTrack* track,
399  LxSimplePoint points[LXSTATIONS][LXLAYERS],
400  Int_t i,
401  Int_t j) {
402  for (; i >= 0; --i) {
403  for (; j >= 0; --j) {
404  if (!track->muchPoints[i][j].empty())
405  points[i][j] = track->muchPoints[i][j].front();
406  else {
407  if (0 != track->parent)
408  return GetPoints(track->parent, points, i, j);
409  else
410  return false;
411  }
412  }
413 
414  j = LXLAYERS - 1;
415  }
416 
417  return true;
418 }
419 
420 static bool GetPoints2(LxSimpleTrack* track,
421  list<LxSimplePoint> points[LXSTATIONS][LXLAYERS],
422  Int_t i,
423  Int_t j) {
424  for (; i >= 0; --i) {
425  for (; j >= 0; --j) {
426  if (!track->muchPoints[i][j].empty())
427  points[i][j].insert(points[i][j].end(),
428  track->muchPoints[i][j].begin(),
429  track->muchPoints[i][j].end());
430  else {
431  if (0 != track->parent)
432  return GetPoints2(track->parent, points, i, j);
433  else
434  return false;
435  }
436  }
437 
438  j = LXLAYERS - 1;
439  }
440 
441  return true;
442 }
443 
444 void LxTrackAnaSegments::BuildStatistics() {
445  Double_t cutCoeff = 4.0;
446  bool readHistoResult = true;
447  Double_t xDispLeft[LXSTATIONS];
448  Double_t yDispLeft[LXSTATIONS];
449  Double_t xDispRight[LXSTATIONS];
450  Double_t yDispRight[LXSTATIONS];
451  Double_t xDispRL[LXSTATIONS];
452  Double_t yDispRL[LXSTATIONS];
453 
454  for (int i = 0; i < LXSTATIONS; ++i) {
455  readHistoResult &= GetHistoRMS("muchInStationXDispLeft", i, xDispLeft[i]);
456  readHistoResult &= GetHistoRMS("muchInStationYDispLeft", i, yDispLeft[i]);
457  readHistoResult &= GetHistoRMS("muchInStationXDispRight", i, xDispRight[i]);
458  readHistoResult &= GetHistoRMS("muchInStationYDispRight", i, yDispRight[i]);
459  readHistoResult &= GetHistoRMS("muchInStationXDispRL", i, xDispRL[i]);
460  readHistoResult &= GetHistoRMS("muchInStationYDispRL", i, yDispRL[i]);
461  }
462 
463  if (!readHistoResult) return;
464 
465  vector<LxSimpleTrack*>& tracks = owner.allTracks;
466 
467  static Int_t triggerEvents = 0;
468  bool hasPositive = false;
469  bool hasNegative = false;
470 
471  for (vector<LxSimpleTrack*>::iterator i = tracks.begin(); i != tracks.end();
472  ++i) {
473  LxSimpleTrack* track = *i;
474 
475  if (useBgr
476  || (0 > track->motherId
477  && (13 == track->pdgCode || -13 == track->pdgCode)))
478  StatForTrack(track);
479 
480  list<LxSimplePoint> points[LXSTATIONS][LXLAYERS];
481 
482  if (!GetPoints2(track, points, stationsInAlgo - 1, LXLAYERS - 1)) continue;
483 
484  bool checkSign = true;
485  bool aOkS[6] = {false, false, false, false, false, false};
486 
487  for (Int_t j = 0; j < stationsInAlgo; ++j) {
488  bool angleOk = false;
489 
490  for (list<LxSimplePoint>::iterator k = points[j][1].begin();
491  k != points[j][1].end();
492  ++k) {
493  LxSimplePoint p1 = *k;
494  Double_t txEst = p1.x / p1.z;
495  Double_t tyEst = p1.y / p1.z;
496 
497  for (list<LxSimplePoint>::iterator l = points[j][2].begin();
498  l != points[j][2].end();
499  ++l) {
500  LxSimplePoint p2 = *l;
501  Double_t deltaZ = p2.z - p1.z;
502  Double_t rX = p1.x + txEst * deltaZ;
503  Double_t rY = p1.y + tyEst * deltaZ;
504 
505  if (abs(rX - p2.x) <= cutCoeff * xDispRight[j]
506  && abs(rY - p2.y) <= cutCoeff * yDispRight[j]) {
507  angleOk = true;
508  break;
509  }
510  } // l
511 
512  if (angleOk) break;
513  } // k
514 
515  if (angleOk) aOkS[j] = true;
516  } // j
517 
518  int okSts = 0;
519 
520  for (int j = 0; j < 6; ++j) {
521  if (aOkS[j]) ++okSts;
522  } // j
523 
524  if (okSts < stationsInAlgo) checkSign = false;
525 
526  if (!checkSign) continue;
527 
528  for (list<LxSimplePoint>::iterator j = points[0][1].begin();
529  j != points[0][1].end();
530  ++j) {
531  LxSimplePoint p01 = *j;
532  Double_t txEst = p01.x / p01.z;
533 
534  for (list<LxSimplePoint>::iterator k = points[1][1].begin();
535  k != points[1][1].end();
536  ++k) {
537  LxSimplePoint p11 = *k;
538  Double_t tx2 = (p11.x - p01.x) / (p11.z - p01.z);
539  Double_t sign2 = tx2 - txEst;
540 
541  for (list<LxSimplePoint>::iterator l = points[0][0].begin();
542  l != points[0][0].end();
543  ++l) {
544  LxSimplePoint p00 = *l;
545 
546  for (list<LxSimplePoint>::iterator m = points[0][2].begin();
547  m != points[0][2].end();
548  ++m) {
549  LxSimplePoint p02 = *m;
550  Double_t tx1 = (p02.x - p00.x) / (p02.z - p00.z);
551  Double_t sign1 = tx1 - txEst;
552 
553  if (sign1 > 0 && sign2 > 0)
554  hasPositive = true;
555  else if (sign1 < 0 && sign2 < 0)
556  hasNegative = true;
557  } // m
558  } // l
559  } // k
560  } // j
561  } // i
562 
563  if (hasPositive && hasNegative) ++triggerEvents;
564 
565  cout << "LxTrackAnaSegments::BuildStatistics(): triggered: " << triggerEvents
566  << " times" << endl;
567 }
568 
569 struct LxSimpleSegment {
570  LxSimplePoint source;
571  LxSimplePoint end;
572  Double_t tx;
573  Double_t ty;
574 
575  LxSimpleSegment() : tx(0), ty(0) {}
576  LxSimpleSegment(LxSimplePoint s, LxSimplePoint e)
577  : source(s)
578  , end(e)
579  , tx((e.x - s.x) / (e.z - s.z))
580  , ty((e.y - s.y) / (e.z - s.z)) {}
581 };
582 
583 void LxTrackAnaSegments::StatForTrack(LxSimpleTrack* track) {
584  for (Int_t i = 0; i < stationsInAlgo; ++i) {
585  for (Int_t j = 0; j < LXLAYERS; ++j) {
586  if (track->muchPoints[i][j].empty()) return;
587  }
588  }
589 
590  LxSimplePoint p1;
591  LxSimplePoint p2;
592  LxSimplePoint p3;
593  Double_t deltaZ;
594  Double_t deltaZ2;
595  Double_t tx;
596  Double_t tx2;
597  Double_t ty;
598  Double_t ty2;
599  Double_t stTx;
600  Double_t stTy;
601  Double_t stTxP;
602  Double_t stTyP;
603 
604  for (Int_t i = 0; i < LXSTATIONS; ++i) {
605  if (track->muchPoints[i][0].empty() || track->muchPoints[i][1].empty()
606  || track->muchPoints[i][2].empty())
607  continue;
608 
609  p1 = track->muchPoints[i][1].front();
610  Double_t txEst = p1.x / p1.z;
611  Double_t tyEst = p1.y / p1.z;
612 
613  p2 = track->muchPoints[i][0].front();
614  deltaZ = p2.z - p1.z;
615  Double_t xEst = p1.x + txEst * deltaZ;
616  Double_t yEst = p1.y + tyEst * deltaZ;
617  muchInStationXDispLeft[i]->Fill(p2.x - xEst);
618  muchInStationYDispLeft[i]->Fill(p2.y - yEst);
619  tx = (p2.x - p1.x) / deltaZ;
620  ty = (p2.y - p1.y) / deltaZ;
621 
622  p2 = track->muchPoints[i][2].front();
623  deltaZ = p2.z - p1.z;
624  xEst = p1.x + txEst * deltaZ;
625  yEst = p1.y + tyEst * deltaZ;
626  muchInStationXDispRight[i]->Fill(p2.x - xEst);
627  muchInStationYDispRight[i]->Fill(p2.y - yEst);
628 
629  tx2 = (p2.x - p1.x) / deltaZ;
630  ty2 = (p2.y - p1.y) / deltaZ;
631 
632  muchInStationTxBreak[i]->Fill(tx2 - tx);
633  muchInStationTyBreak[i]->Fill(ty2 - ty);
634 
635  stTxP = stTx;
636  stTyP = stTy;
637  p1 = track->muchPoints[i][0].front();
638  deltaZ = p2.z - p1.z;
639  stTx = (p2.x - p1.x) / deltaZ;
640  stTy = (p2.y - p1.y) / deltaZ;
641  muchInStationXDispRL[i]->Fill(p1.x - p2.x + tx2 * deltaZ);
642  muchInStationYDispRL[i]->Fill(p1.y - p2.y + ty2 * deltaZ);
643 
644  if (i > 0) {
645  p1 = track->muchPoints[i - 1][LXMIDDLE].front();
646  p2 = track->muchPoints[i][LXMIDDLE].front();
647  deltaZ = p2.z - p1.z;
648  tx = (p2.x - p1.x) / deltaZ;
649  muchLongSegmentTxHisto[i - 1]->Fill(tx - p2.x / p2.z);
650  ty = (p2.y - p1.y) / deltaZ;
651  muchLongSegmentTyHisto[i - 1]->Fill(ty - p2.y / p2.z);
652 
653  muchOutStationTxBreakRight[i - 1]->Fill(tx - stTx);
654  muchOutStationTyBreakRight[i - 1]->Fill(ty - stTy);
655 
656  muchOutStationTxBreakLeft[i - 1]->Fill(tx - stTxP);
657  muchOutStationTyBreakLeft[i - 1]->Fill(ty - stTyP);
658 
659  muchOutStationXDispByTriplet[i - 1]->Fill(p1.x - p2.x + stTx * deltaZ);
660  muchOutStationYDispByTriplet[i - 1]->Fill(p1.y - p2.y + stTy * deltaZ);
661 
662  Double_t txVertex = p2.x / p2.z;
663  Double_t tyVertex = p2.y / p2.z;
664  //Double_t scatCoeff = sqrt(1 + txVertex * txVertex + tyVertex * tyVertex);
665 
666  muchOutStationXDispByVertex[i - 1]->Fill(
667  p1.x - p2.x + txVertex * deltaZ /* / scatCoeff*/);
668  muchOutStationYDispByVertex[i - 1]->Fill(
669  p1.y - p2.y + tyVertex * deltaZ /* / scatCoeff*/);
670 
671  muchTripletTxBreak[i - 1]->Fill(stTxP - stTx);
672  muchTripletTyBreak[i - 1]->Fill(stTyP - stTy);
673 
674  // Rather complex part for implementation: calculate the dispersion characteristics for segment clusters.
675  Double_t maxXdisp = 0;
676  Double_t maxYdisp = 0;
677  Double_t maxTxdisp = 0;
678  Double_t maxTydisp = 0;
679 
680  for (list<LxSimplePoint>::iterator l0 =
681  track->muchPoints[i - 1][0].begin();
682  l0 != track->muchPoints[i - 1][0].end();
683  ++l0) {
684  for (list<LxSimplePoint>::iterator l1 =
685  track->muchPoints[i - 1][1].begin();
686  l1 != track->muchPoints[i - 1][1].end();
687  ++l1) {
688  for (list<LxSimplePoint>::iterator l2 =
689  track->muchPoints[i - 1][2].begin();
690  l2 != track->muchPoints[i - 1][2].end();
691  ++l2) {
692  for (list<LxSimplePoint>::iterator r0 =
693  track->muchPoints[i][0].begin();
694  r0 != track->muchPoints[i][0].end();
695  ++r0) {
696  for (list<LxSimplePoint>::iterator r1 =
697  track->muchPoints[i][1].begin();
698  r1 != track->muchPoints[i][1].end();
699  ++r1) {
700  for (list<LxSimplePoint>::iterator r2 =
701  track->muchPoints[i][2].begin();
702  r2 != track->muchPoints[i][2].end();
703  ++r2) {
704  LxSimplePoint lPoints[LXLAYERS] = {*l0, *l1, *l2};
705  LxSimplePoint rPoints[LXLAYERS] = {*r0, *r1, *r2};
706  LxSimpleSegment segments[LXLAYERS * LXLAYERS] = {};
707 
708  for (Int_t j = 0; j < LXLAYERS; ++j) {
709  for (Int_t k = 0; k < LXLAYERS; ++k)
710  segments[j * LXLAYERS + k] =
711  LxSimpleSegment(rPoints[j], lPoints[k]);
712  }
713 
714  for (Int_t j = 0; j < LXLAYERS * LXLAYERS - 1; ++j) {
715  LxSimpleSegment s1 = segments[j];
716 
717  for (Int_t k = j + 1; k < LXLAYERS * LXLAYERS; ++k) {
718  LxSimpleSegment s2 = segments[k];
719  Double_t diffZ = s1.source.z - s2.source.z;
720  Double_t dtx = abs(s2.tx - s1.tx);
721  Double_t dty = abs(s2.ty - s1.ty);
722  Double_t dx =
723  abs(s2.source.x + s2.tx * diffZ - s1.source.x);
724  Double_t dy =
725  abs(s2.source.y + s2.ty * diffZ - s1.source.y);
726 
727  if (maxXdisp < dx) maxXdisp = dx;
728 
729  if (maxYdisp < dy) maxYdisp = dy;
730 
731  if (maxTxdisp < dtx) maxTxdisp = dtx;
732 
733  if (maxTydisp < dty) maxTydisp = dty;
734  }
735  }
736  }
737  }
738  }
739  }
740  }
741  } // for (list<LxSimplePoint>::iterator l0 = track->muchPoints[i - 1][0].begin(); l0 != track->muchPoints[i - 1][0].end(); ++l0)
742 
743  muchClusterXDispHisto[i - 1]->Fill(maxXdisp);
744  muchClusterYDispHisto[i - 1]->Fill(maxYdisp);
745  muchClusterTxDispHisto[i - 1]->Fill(maxTxdisp);
746  muchClusterTyDispHisto[i - 1]->Fill(maxTydisp);
747 
748  if (i < LXSTATIONS - 1) {
749  p1 = track->muchPoints[i - 1][LXMIDDLE].front();
750  p2 = track->muchPoints[i][LXMIDDLE].front();
751  deltaZ = p2.z - p1.z;
752  tx = (p2.x - p1.x) / deltaZ;
753  ty = (p2.y - p1.y) / deltaZ;
754  p3 = track->muchPoints[i + 1][LXMIDDLE].front();
755  deltaZ2 = p3.z - p2.z;
756  tx2 = (p3.x - p2.x) / deltaZ2;
757  Double_t slopeCoeff = sqrt(1 + tx2 * tx2 + ty2 * ty2);
758  muchSegmentTxBreakHisto[i - 1]->Fill((tx2 - tx) / slopeCoeff);
759  ty2 = (p3.y - p2.y) / deltaZ2;
760  muchSegmentTyBreakHisto[i - 1]->Fill((ty2 - ty) / slopeCoeff);
761  }
762 
763  for (Int_t j = 0; j < LXLAYERS * LXLAYERS; ++j) {
764  p1 = track->muchPoints[i - 1][j % LXLAYERS].front();
765  p2 = track->muchPoints[i][j / LXLAYERS].front();
766  deltaZ = p2.z - p1.z;
767  tx = (p2.x - p1.x) / deltaZ;
768  ty = (p2.y - p1.y) / deltaZ;
769 
770  for (Int_t k = j + 1; k < LXLAYERS * LXLAYERS; ++k) {
771  p1 = track->muchPoints[i - 1][k % LXLAYERS].front();
772  p2 = track->muchPoints[i][k / LXLAYERS].front();
773  deltaZ = p2.z - p1.z;
774  tx2 = (p2.x - p1.x) / deltaZ;
775  muchStationTxDispHisto[i - 1]->Fill(tx2 - tx);
776  ty2 = (p2.y - p1.y) / deltaZ;
777  muchStationTyDispHisto[i - 1]->Fill(ty2 - ty);
778  }
779  }
780  }
781 
782  if (i < LXSTATIONS - 1) {
783  p1 = track->muchPoints[i + 1][LXMIDDLE].front();
784  p2 = track->muchPoints[i][LXMIDDLE].front();
785  deltaZ = p2.z - p1.z;
786  Double_t deltaX = p2.x - p1.x - p1.tx * deltaZ;
787  Double_t deltaY = p2.y - p1.y - p1.ty * deltaZ;
788  Double_t deltaTx = p2.tx - p1.tx;
789  Double_t deltaTy = p2.ty - p1.ty;
790  muchXTxCovHisto[i]->Fill(deltaX, deltaTx);
791  muchYTyCovHisto[i]->Fill(deltaY, deltaTy);
792  }
793  }
794 }
muchInStationXDispRL
static TH1F * muchInStationXDispRL[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:20
LxSimpleSegment::ty
scaltype ty
Definition: Simple/LxTrackAnaSegments.cxx:369
LxSimpleTrack::motherId
Int_t motherId
Definition: Simple/LxTrackAna.h:27
GetPoints
static bool GetPoints(LxSimpleTrack *track, LxSimplePoint points[LXSTATIONS][LXLAYERS], Int_t i, Int_t j)
Definition: riplet/LxTrackAnaSegments.cxx:398
muchOutStationTyBreakRight
static TH1F * muchOutStationTyBreakRight[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:29
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
muchOutStationYDispByVertex
static TH1F * muchOutStationYDispByVertex[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:38
muchOutStationYDispByTriplet
static TH1F * muchOutStationYDispByTriplet[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:35
muchTripletTyBreak
static TH1F * muchTripletTyBreak[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:32
SaveHisto
static void SaveHisto(TH1 *histo)
Definition: riplet/LxTrackAnaSegments.cxx:329
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
LxTrackAnaSegments::owner
LxTrackAna & owner
Definition: Simple/LxTrackAnaSegments.h:21
muchInStationTyBreak
static TH1F * muchInStationTyBreak[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:24
LxTrackAnaSegments::SetParticleType
void SetParticleType(TString v)
Definition: Simple/LxTrackAnaSegments.cxx:56
LxSimpleSegment::source
LxSimplePoint source
Definition: Simple/LxTrackAnaSegments.cxx:366
muchStationTyDispHisto
static TH1F * muchStationTyDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:47
saveHistos
bool saveHistos
Definition: riplet/LxTrackAnaSegments.cxx:11
LxSimpleSegment::tx
Double_t tx
Definition: riplet/LxTrackAnaSegments.cxx:572
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
muchSegmentTyBreakHisto
static TH1F * muchSegmentTyBreakHisto[LXSTATIONS - 2]
Definition: riplet/LxTrackAnaSegments.cxx:44
muchStationTxDispHisto
static TH1F * muchStationTxDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:46
muchInStationXDispRight
static TH1F * muchInStationXDispRight[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:16
LXLAYERS
#define LXLAYERS
Definition: Simple/LxSettings.h:8
LxSimplePoint::y
scaltype y
Definition: Simple/LxTrackAna.h:16
muchInStationXDispLeft
static TH1F * muchInStationXDispLeft[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:15
muchYTyCovHisto
static TH2F * muchYTyCovHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:50
LXMIDDLE
#define LXMIDDLE
Definition: Simple/LxSettings.h:10
LxSimpleTrack
Definition: Simple/LxTrackAna.h:25
LxTrackAnaTriplet
Definition: riplet/LxTrackAna.h:68
LxSimpleSegment::ty
Double_t ty
Definition: riplet/LxTrackAnaSegments.cxx:573
muchOutStationTxBreakRight
static TH1F * muchOutStationTxBreakRight[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:27
LxTrackAna.h
LxSimpleSegment::tx
scaltype tx
Definition: Simple/LxTrackAnaSegments.cxx:368
LxSimplePoint::z
scaltype z
Definition: Simple/LxTrackAna.h:17
tracks
TClonesArray * tracks
Definition: Analyze_matching.h:17
LxTrackAnaSegments::BuildStatistics
void BuildStatistics()
Definition: Simple/LxTrackAnaSegments.cxx:353
muchInStationYDispLeft
static TH1F * muchInStationYDispLeft[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:17
LxSimpleSegment
Definition: Simple/LxTrackAnaSegments.cxx:365
LxTrackAna::allTracks
std::vector< LxSimpleTrack * > allTracks
Definition: Simple/LxTrackAna.h:118
LxTrackAnaSegments::StatForTrack
void StatForTrack(LxSimpleTrack *track)
Definition: Simple/LxTrackAnaSegments.cxx:379
muchOutStationXDispByTriplet
static TH1F * muchOutStationXDispByTriplet[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:34
LxSimplePoint::x
scaltype x
Definition: Simple/LxTrackAna.h:15
muchClusterYDispHisto
static TH1F * muchClusterYDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:53
LxTrackAnaSegments::stationsInAlgo
Int_t stationsInAlgo
Definition: Simple/LxTrackAnaSegments.h:22
cutCoeff
Double_t cutCoeff
Definition: riplet/Lx.cxx:60
muchLongSegmentTyHisto
static TH1F * muchLongSegmentTyHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:41
muchInStationYDispRight
static TH1F * muchInStationYDispRight[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:18
muchOutStationXDispByVertex
static TH1F * muchOutStationXDispByVertex[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:37
LxTrackAnaSegments::useBgr
bool useBgr
Definition: riplet/LxTrackAnaSegments.h:29
muchClusterTxDispHisto
static TH1F * muchClusterTxDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:54
LxTrackAnaSegments::Init
void Init()
Definition: Simple/LxTrackAnaSegments.cxx:62
LxTrackAnaSegments.h
muchSegmentTxBreakHisto
static TH1F * muchSegmentTxBreakHisto[LXSTATIONS - 2]
Definition: riplet/LxTrackAnaSegments.cxx:43
GetHistoRMS
static bool GetHistoRMS(const char *histoNameBase, Int_t histoNumber, Double_t &retVal)
Definition: riplet/LxTrackAnaSegments.cxx:58
points
TClonesArray * points
Definition: Analyze_matching.h:18
LxSimpleTrack::parent
LxSimpleTrack * parent
Definition: riplet/LxTrackAna.h:64
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
LxSimplePoint::ty
scaltype ty
Definition: Simple/LxTrackAna.h:19
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
muchInStationYDispRL
static TH1F * muchInStationYDispRL[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:21
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
muchLongSegmentTxHisto
static TH1F * muchLongSegmentTxHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:40
LxSimpleSegment::end
LxSimplePoint end
Definition: Simple/LxTrackAnaSegments.cxx:367
LxSimpleTrack::muchPoints
std::list< LxSimplePoint > muchPoints[LXSTATIONS][LXLAYERS]
Definition: Simple/LxTrackAna.h:55
LxTrackAnaSegments::LxTrackAnaSegments
LxTrackAnaSegments(LxTrackAna &o)
Definition: Simple/LxTrackAnaSegments.cxx:51
LxSimplePoint::tx
scaltype tx
Definition: Simple/LxTrackAna.h:18
GetPoints2
static bool GetPoints2(LxSimpleTrack *track, list< LxSimplePoint > points[LXSTATIONS][LXLAYERS], Int_t i, Int_t j)
Definition: riplet/LxTrackAnaSegments.cxx:420
muchOutStationTxBreakLeft
static TH1F * muchOutStationTxBreakLeft[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:26
muchInStationTxBreak
static TH1F * muchInStationTxBreak[LXSTATIONS]
Definition: riplet/LxTrackAnaSegments.cxx:23
muchOutStationTyBreakLeft
static TH1F * muchOutStationTyBreakLeft[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:28
LxSimpleSegment::LxSimpleSegment
LxSimpleSegment()
Definition: riplet/LxTrackAnaSegments.cxx:575
LxTrackAnaSegments::Finish
void Finish()
Definition: Simple/LxTrackAnaSegments.cxx:307
LXSTATIONS
#define LXSTATIONS
Definition: Simple/LxSettings.h:9
LxSimpleTrack::pdgCode
Int_t pdgCode
Definition: Simple/LxTrackAna.h:26
muchTripletTxBreak
static TH1F * muchTripletTxBreak[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:31
particleType
static TString particleType("jpsi")
LxSimpleSegment::LxSimpleSegment
LxSimpleSegment(LxSimplePoint s, LxSimplePoint e)
Definition: riplet/LxTrackAnaSegments.cxx:576
muchClusterXDispHisto
static TH1F * muchClusterXDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:52
muchXTxCovHisto
static TH2F * muchXTxCovHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:49
muchClusterTyDispHisto
static TH1F * muchClusterTyDispHisto[LXSTATIONS - 1]
Definition: riplet/LxTrackAnaSegments.cxx:55
LxSimplePoint
Definition: Simple/LxTrackAna.h:14