1 |
pam-fi |
1.1 |
/* |
2 |
|
|
* TrkDedxHeCut.cpp |
3 |
|
|
* |
4 |
pam-fi |
1.2 |
* Created on: 28/gen/2010 |
5 |
|
|
* Author: Nicola Mori |
6 |
pam-fi |
1.1 |
*/ |
7 |
|
|
|
8 |
|
|
/*! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. */ |
9 |
|
|
|
10 |
|
|
#include "TrkDedxHeCut.h" |
11 |
|
|
|
12 |
|
|
int TrkDedxHeCut::Check(PamLevel2 *event) { |
13 |
|
|
|
14 |
pam-fi |
1.2 |
TrkTrack *track = event->GetTrack(0)->GetTrkTrack(); |
15 |
|
|
float R = track->GetRigidity(); |
16 |
|
|
//float dEdx = track->GetDEDX(); |
17 |
|
|
|
18 |
|
|
// New smart dE/dx computation |
19 |
|
|
// This procedure removes high releases and computes dE/dx |
20 |
|
|
float totReleaseX = 0., maxReleaseX = 0.; |
21 |
|
|
float totReleaseY = 0., maxReleaseY = 0.; |
22 |
|
|
int nX = 0, nY = 0; |
23 |
|
|
float dEdxView = 0; |
24 |
|
|
/* |
25 |
|
|
ToFLevel2 *tofL2 = event->GetToFLevel2(); |
26 |
|
|
int trkSeqNo = 0; |
27 |
|
|
for (int i = 0; i < event->GetToFLevel2()->ntrk(); i++) { |
28 |
|
|
if (tofL2->GetToFTrkVar(i)->trkseqno == track->GetSeqNo()) { |
29 |
|
|
trkSeqNo = i; |
30 |
|
|
break; |
31 |
|
|
} |
32 |
|
|
} |
33 |
|
|
if ((tofL2->GetdEdx(trkSeqNo, 0, 100) + tofL2->GetdEdx(trkSeqNo, 1, 100) + tofL2->GetdEdx(trkSeqNo, 2, 100) |
34 |
|
|
+ tofL2->GetdEdx(trkSeqNo, 3, 100)) / 4. < 2) { |
35 |
|
|
cout << "Check" << endl; |
36 |
|
|
for (int ip = 0; ip < 6; ip++) { |
37 |
|
|
cout << "TOF" << ip << ": " << tofL2->GetdEdx(trkSeqNo, ip, 100) << endl; |
38 |
|
|
dEdxView = track->GetDEDX(ip, 0); |
39 |
|
|
if (dEdxView > 0 && track->XGood(ip)) |
40 |
|
|
cout << "X" << ip << ": " << dEdxView << endl; |
41 |
|
|
dEdxView = track->GetDEDX(ip, 1); |
42 |
|
|
if (dEdxView > 0 && track->YGood(ip)) |
43 |
|
|
cout << "Y" << ip << ": " << dEdxView << endl; |
44 |
|
|
} |
45 |
|
|
} |
46 |
|
|
*/ |
47 |
|
|
for (int ip = 0; ip < 6; ip++) { |
48 |
|
|
// X view |
49 |
|
|
dEdxView = track->GetDEDX(ip, 0); |
50 |
|
|
if (dEdxView > 0 && track->XGood(ip)) { |
51 |
|
|
totReleaseX += dEdxView; |
52 |
|
|
nX++; |
53 |
|
|
if (dEdxView > maxReleaseX) |
54 |
|
|
maxReleaseX = dEdxView; |
55 |
|
|
} |
56 |
|
|
// Y view |
57 |
|
|
dEdxView = track->GetDEDX(ip, 1); |
58 |
|
|
if (dEdxView > 0 && track->YGood(ip)) { |
59 |
|
|
totReleaseY += dEdxView; |
60 |
|
|
nY++; |
61 |
|
|
if (dEdxView > maxReleaseY) |
62 |
|
|
maxReleaseY = dEdxView; |
63 |
|
|
} |
64 |
|
|
} |
65 |
pam-fi |
1.1 |
|
66 |
pam-fi |
1.2 |
// Discard highest release, eventually |
67 |
|
|
if (maxReleaseX > -9. + 4. * totReleaseX / (float) nX) { |
68 |
|
|
totReleaseX -= maxReleaseX; |
69 |
|
|
nX--; |
70 |
|
|
} |
71 |
|
|
if (maxReleaseY > -8. + 4. * totReleaseY / (float) nY) { |
72 |
|
|
totReleaseY -= maxReleaseY; |
73 |
|
|
nY--; |
74 |
pam-fi |
1.1 |
} |
75 |
|
|
|
76 |
pam-fi |
1.2 |
// Compute dE/dx |
77 |
|
|
//float dEdx = (totReleaseX + totReleaseY) / (float) (nX + nY); |
78 |
|
|
float dEdxX = totReleaseX / nX; |
79 |
|
|
float dEdxY = totReleaseY / nY; |
80 |
|
|
|
81 |
|
|
// Analyze the event |
82 |
|
|
// X |
83 |
|
|
float denLow = pow(R, 1.8); // The power of the denominator is the same for X and Y (with current calibration) |
84 |
|
|
float denHigh = pow(R, 1.5); |
85 |
|
|
if (dEdxX < 3.7 + 4.6 / denLow) |
86 |
|
|
return 0; |
87 |
|
|
if (dEdxX > 8.9 + 17. / denHigh) |
88 |
|
|
return 0; |
89 |
|
|
// Y |
90 |
|
|
if (dEdxY < 3.3 + 4.9 / denLow) |
91 |
|
|
return 0; |
92 |
|
|
if (dEdxY > 8.0 + 17. / denHigh) |
93 |
|
|
return 0; |
94 |
|
|
//X+Y |
95 |
|
|
/*if (dEdx < 3.7 + 4.8 / (R * R)) |
96 |
|
|
return 0; |
97 |
|
|
if (dEdx > 8.9 + 17. / (R * R)) |
98 |
|
|
return 0;*/ |
99 |
pam-fi |
1.1 |
|
100 |
pam-fi |
1.2 |
return CUTOK; |
101 |
pam-fi |
1.1 |
} |