/[PAMELA software]/PamCut/TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp
ViewVC logotype

Annotation of /PamCut/TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Mon May 3 14:36:27 2010 UTC (14 years, 7 months ago) by pam-fi
Branch: MAIN
CVS Tags: Root_V8, MergedToHEAD_1, nuclei_reproc, MergedFromV8_1, BeforeMergingFromV8_1, V9
Branch point for: V8
Changes since 1.1: +85 -51 lines
Cut changed: new function, independent check on X and Y views, improved proton rejection with maximum-release cluster elimination.

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 }

  ViewVC Help
Powered by ViewVC 1.1.23