TrkCuts/TrkNucleiZCut/TrkNucleiZCut.cpp

/*
 * TrkNucleiZCut.cpp
 *
 *  Created on: 10-mag-2009
 *      Author: Nicola Mori
 */

/*! @file TrkNucleiZCut.cpp The TrkNucleiZCut class implementation file */
#ifndef NO_TRKNUCLEI

#include "TrkNucleiZCut.h"

const float TrkNucleiZCut::_meanRigMean[] = { .992, 1.99, 3., 4., 5., 6. }; // the charge distribution peak for Z=1..6
const float TrkNucleiZCut::_sigmaRigMean[] = { 0.06, 0.1, 0.17, 0.34, 0.44, 0.59 }; // the charge-width for Z=1..6

const float TrkNucleiZCut::_meanBetaMean[] = { 1., 2., 3., 4., 5., 6. }; // the charge distribution peak for Z=1..6
const float TrkNucleiZCut::_sigmaBetaMean[] = { 0.09, 0.1, 0.17, 0.34, 0.44, 0.59 }; // the charge-width for Z=1..6


int TrkNucleiZCut::Check(PamLevel2 *event) {

  if (_lowZ < 1 || _lowZ > 6) {
    return TRKNUCLEIZ_ILLEGALLOWZ;
  }
  if ((_highZ < (int) _lowZ && _highZ != -1) || _highZ > 6) {
    return TRKNUCLEIZ_ILLEGALLOWZ;
  }

  /* Check if we must perform a range selection */
  bool range = false;
  if (_highZ != -1)
    range = true;

  // ******** Check if TrkNuclei has already been initialized ********
  if (_trkNuclei == NULL)
    _trkNuclei = new TrkNuclei();
  _trkNuclei->Set(event->GetTrack(0)->GetTrkTrack());

  // ******** First of all, check if the minimum number of planes is matched. ********
  // NOTE: in current implementation a plane is considered valid if it is tagged as good in Tracker level 2
  //       routines. This should correspond to a valid charge.
  // TODO Check the above statement.

  // Prepare the valid charges counter
  unsigned int nValid = 0;
  for (unsigned int ip = 0; ip < 6; ip++) {
    if (_trkNuclei->GetTrkTrack()->XGood(ip) && _trkNuclei->GetTrkTrack()->XGood(ip))
      nValid++;
  }
  if (nValid < _minLayers)
    return TRKNUCLEIZ_TOOFEWLAYERS;

  // ******** Check the mean charge value ********
  if (_howToCheck == TRKNUCLEIZ_CHECKMEAN) {

    Float_t charge = 0., meanLowZ = 0., sigmaLowZ = 0., meanHighZ = 0., sigmaHighZ = 0.;
    ;
    switch (_method) {
    case TRKNUCLEIZ_RIG:
      charge = _trkNuclei->GetZ_Rigidity();
      meanLowZ = _meanRigMean[_lowZ - 1];
      sigmaLowZ = _sigmaRigMean[_lowZ - 1];
      if (range) {
        meanHighZ = _meanRigMean[_highZ - 1];
        sigmaHighZ = _sigmaRigMean[_highZ - 1];
      }
      break;
    case TRKNUCLEIZ_BETA:
      charge = _trkNuclei->GetZ_Beta(event->GetTrack(0)->GetToFTrack()->beta[12]);
      meanLowZ = _meanBetaMean[_lowZ - 1];
      sigmaLowZ = _sigmaBetaMean[_lowZ - 1];
      if (range) {
        meanHighZ = _meanBetaMean[_highZ - 1];
        sigmaHighZ = _sigmaBetaMean[_highZ - 1];
      }

      break;
    }

    if (!range) {
      if ((charge > meanLowZ + _upperLimit * sigmaLowZ) || (charge < meanLowZ - _lowerLimit * sigmaLowZ)) {
        return TRKNUCLEIZ_OUTOFBOUNDS;
      }
    }
    else {
      if ((charge > meanHighZ + _upperLimit * sigmaHighZ) || (charge < meanLowZ - _lowerLimit * sigmaLowZ)) {
        return TRKNUCLEIZ_OUTOFBOUNDS;
      }
    }
  }
  // ******** Check the charge value from each plane ********
  // NB. In current implementation only 2 check methods (mean and single values) are used, so an "else" statement
  //     is sufficient. Implement an "if" cascade or a "switch" if another check method is necessary.

  else {
    Float_t charge[6] = { 0., 0., 0., 0., 0., 0. };
    nValid = 0;
    switch (_method) {
    case TRKNUCLEIZ_RIG:
      for (int ip = 0; ip < 6; ip++) {
        charge[ip] = _trkNuclei->GetZ_Rigidity(ip);
        if (charge[ip] > 0.) {
          nValid++;
        }
      }
      break;
    case TRKNUCLEIZ_BETA:
      for (int ip = 0; ip < 6; ip++) {
        charge[ip] = _trkNuclei->GetZ_Beta(ip, event->GetTrack(0)->GetToFTrack()->beta[12]);
        if (charge[ip] > 0.) {
          nValid++;
        }
      }
      break;
    }

    for (int ip = 0; ip < 6; ip++) {
      if (charge[ip] > 0.) { // Check only good layers
        // Currently, no calibration for single layers is available. Z is used as a mean value, while for sigma
        // we use the standard values use for mean dE/dx Vs. beta.
        if (!range) {
          if ((charge[ip] > _lowZ + _upperLimit * _sigmaBetaMean[_lowZ - 1]) || (charge[ip] < _lowZ - _lowerLimit
              * _sigmaBetaMean[_lowZ - 1])) {
            return TRKNUCLEIZ_OUTOFBOUNDS;
          }
        }
        else {
          if ((charge[ip] > _highZ + _upperLimit * _sigmaBetaMean[_highZ - 1]) || (charge[ip] < _lowZ - _lowerLimit
              * _sigmaBetaMean[_lowZ - 1])) {
            return TRKNUCLEIZ_OUTOFBOUNDS;
          }
        }
      }
    }
  }

  return CUTOK;

}

#endif /* NO_TRKNUCLEI */
1	pam-fi	1.1	/*
2			* TrkNucleiZCut.cpp
3			*
4			* Created on: 10-mag-2009
5			* Author: Nicola Mori
6			*/
7
8			/! @file TrkNucleiZCut.cpp The TrkNucleiZCut class implementation file /
9			#ifndef NO_TRKNUCLEI
10
11			#include "TrkNucleiZCut.h"
12
13	pam-fi	1.4	const float TrkNucleiZCut::_meanRigMean[] = { .992, 1.99, 3., 4., 5., 6. }; // the charge distribution peak for Z=1..6
14			const float TrkNucleiZCut::_sigmaRigMean[] = { 0.06, 0.1, 0.17, 0.34, 0.44, 0.59 }; // the charge-width for Z=1..6
15	pam-fi	1.2
16	pam-fi	1.4	const float TrkNucleiZCut::_meanBetaMean[] = { 1., 2., 3., 4., 5., 6. }; // the charge distribution peak for Z=1..6
17			const float TrkNucleiZCut::_sigmaBetaMean[] = { 0.09, 0.1, 0.17, 0.34, 0.44, 0.59 }; // the charge-width for Z=1..6
18	pam-fi	1.1
19
20			int TrkNucleiZCut::Check(PamLevel2 *event) {
21
22	pam-fi	1.3	if (_lowZ < 1 \|\| _lowZ > 6) {
23			return TRKNUCLEIZ_ILLEGALLOWZ;
24			}
25	pam-fi	1.4	if ((_highZ < (int) _lowZ && _highZ != -1) \|\| _highZ > 6) {
26	pam-fi	1.3	return TRKNUCLEIZ_ILLEGALLOWZ;
27			}
28
29			/* Check if we must perform a range selection */
30			bool range = false;
31	pam-fi	1.4	if (_highZ != -1)
32			range = true;
33	pam-fi	1.1
34			// ****** Check if TrkNuclei has already been initialized ******
35			if (_trkNuclei == NULL)
36	pam-fi	1.4	_trkNuclei = new TrkNuclei();
37	pam-fi	1.1	_trkNuclei->Set(event->GetTrack(0)->GetTrkTrack());
38
39			// ****** First of all, check if the minimum number of planes is matched. ******
40			// NOTE: in current implementation a plane is considered valid if it is tagged as good in Tracker level 2
41			// routines. This should correspond to a valid charge.
42			// TODO Check the above statement.
43
44			// Prepare the valid charges counter
45			unsigned int nValid = 0;
46			for (unsigned int ip = 0; ip < 6; ip++) {
47			if (_trkNuclei->GetTrkTrack()->XGood(ip) && _trkNuclei->GetTrkTrack()->XGood(ip))
48	pam-fi	1.4	nValid++;
49	pam-fi	1.1	}
50			if (nValid < _minLayers)
51	pam-fi	1.4	return TRKNUCLEIZ_TOOFEWLAYERS;
52	pam-fi	1.1
53			// ****** Check the mean charge value ******
54			if (_howToCheck == TRKNUCLEIZ_CHECKMEAN) {
55
56	pam-fi	1.4	Float_t charge = 0., meanLowZ = 0., sigmaLowZ = 0., meanHighZ = 0., sigmaHighZ = 0.;
57			;
58	pam-fi	1.1	switch (_method) {
59	pam-fi	1.4	case TRKNUCLEIZ_RIG:
60	pam-fi	1.1	charge = _trkNuclei->GetZ_Rigidity();
61	pam-fi	1.4	meanLowZ = _meanRigMean[_lowZ - 1];
62			sigmaLowZ = _sigmaRigMean[_lowZ - 1];
63			if (range) {
64			meanHighZ = _meanRigMean[_highZ - 1];
65			sigmaHighZ = _sigmaRigMean[_highZ - 1];
66	pam-fi	1.3	}
67	pam-fi	1.1	break;
68	pam-fi	1.4	case TRKNUCLEIZ_BETA:
69	pam-fi	1.1	charge = _trkNuclei->GetZ_Beta(event->GetTrack(0)->GetToFTrack()->beta[12]);
70	pam-fi	1.4	meanLowZ = _meanBetaMean[_lowZ - 1];
71			sigmaLowZ = _sigmaBetaMean[_lowZ - 1];
72			if (range) {
73			meanHighZ = _meanBetaMean[_highZ - 1];
74			sigmaHighZ = _sigmaBetaMean[_highZ - 1];
75	pam-fi	1.3	}
76
77	pam-fi	1.1	break;
78			}
79
80	pam-fi	1.4	if (!range) {
81	pam-fi	1.3	if ((charge > meanLowZ + _upperLimit * sigmaLowZ) \|\| (charge < meanLowZ - _lowerLimit * sigmaLowZ)) {
82			return TRKNUCLEIZ_OUTOFBOUNDS;
83			}
84			}
85			else {
86			if ((charge > meanHighZ + _upperLimit * sigmaHighZ) \|\| (charge < meanLowZ - _lowerLimit * sigmaLowZ)) {
87			return TRKNUCLEIZ_OUTOFBOUNDS;
88			}
89	pam-fi	1.1	}
90			}
91			// ****** Check the charge value from each plane ******
92			// NB. In current implementation only 2 check methods (mean and single values) are used, so an "else" statement
93			// is sufficient. Implement an "if" cascade or a "switch" if another check method is necessary.
94
95			else {
96	pam-fi	1.4	Float_t charge[6] = { 0., 0., 0., 0., 0., 0. };
97	pam-fi	1.1	nValid = 0;
98			switch (_method) {
99	pam-fi	1.4	case TRKNUCLEIZ_RIG:
100	pam-fi	1.1	for (int ip = 0; ip < 6; ip++) {
101			charge[ip] = _trkNuclei->GetZ_Rigidity(ip);
102			if (charge[ip] > 0.) {
103			nValid++;
104			}
105			}
106			break;
107	pam-fi	1.4	case TRKNUCLEIZ_BETA:
108	pam-fi	1.1	for (int ip = 0; ip < 6; ip++) {
109			charge[ip] = _trkNuclei->GetZ_Beta(ip, event->GetTrack(0)->GetToFTrack()->beta[12]);
110			if (charge[ip] > 0.) {
111			nValid++;
112			}
113			}
114			break;
115			}
116
117			for (int ip = 0; ip < 6; ip++) {
118	pam-fi	1.4	if (charge[ip] > 0.) { // Check only good layers
119			// Currently, no calibration for single layers is available. Z is used as a mean value, while for sigma
120			// we use the standard values use for mean dE/dx Vs. beta.
121			if (!range) {
122			if ((charge[ip] > _lowZ + _upperLimit * _sigmaBetaMean[_lowZ - 1]) \|\| (charge[ip] < _lowZ - _lowerLimit
123			* _sigmaBetaMean[_lowZ - 1])) {
124			return TRKNUCLEIZ_OUTOFBOUNDS;
125			}
126	pam-fi	1.1	}
127	pam-fi	1.4	else {
128			if ((charge[ip] > _highZ + _upperLimit * _sigmaBetaMean[_highZ - 1]) \|\| (charge[ip] < _lowZ - _lowerLimit
129			* _sigmaBetaMean[_lowZ - 1])) {
130			return TRKNUCLEIZ_OUTOFBOUNDS;
131			}
132	pam-fi	1.3	}
133			}
134	pam-fi	1.1	}
135			}
136
137			return CUTOK;
138
139			}
140
141			#endif /* NO_TRKNUCLEI */