TrkCuts/TrkDedxHCut/TrkDedxHCut.cpp

/*
 * TrkDedxHCut.cpp
 *
 *  Created on: 15-mar-2010
 *      Author: S. Ricciarini
 */

/*! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file */

#include "TrkDedxHCut.h"

const Double_t TrkDedxHCut::_xLow[] = { 0.126414, 0.140142, 0.152714, 0.165854, 0.178357, 0.194244, 0.211669, 0.242859,
    0.264646, 0.293385, 0.314258, 0.360565, 0.392911, 0.450807, 0.508423, 0.603735, 0.680897, 0.754838, 0.851312,
    1.00032, 1.36201, 1.37, 1000. }; // _nLow=23 elements

const Double_t TrkDedxHCut::_yLow[] = { 18.9221, 16.0752, 14.2214, 12.6324, 11.2421, 9.71928, 8.1965, 6.93856, 6.14407,
    5.15095, 4.6875, 3.8268, 3.49576, 2.89989, 2.56886, 2.17161, 1.90678, 1.70816, 1.50953, 1.27754, 0.715042, 0, 0 }; // _nLow=23 elements

Int_t TrkDedxHCut::Check(PamLevel2 *event) {

//  TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
    if(event->GetNTracks(_trkAlg)==0)return 0;    
    ExtTrack *trkTrack = event->GetTrack(0,_trkAlg)->GetExtTrack();

  Float_t rigMod = trkTrack->GetRigidity();
  Float_t dedx = trkTrack->GetDEDX();

  if (rigMod < 0.126414)
    return 0; // a proton below 126.414 MV is rejected

  // lower limit taken from N. De Simone (interpolation of _xLow)
  
  for (Int_t i = 0; i < _nLow - 1; i++) {
    if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
      Double_t slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
      Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);

      if (dedx < y)
        return 0;
    }
  }

  // higher limit taken from N. Mori lower limit for Helium

  // Compute dE/dx for X and Y separately
  Float_t totReleaseX = 0.;
  Float_t totReleaseY = 0.;
  Int_t nX = 0, nY = 0;
  Float_t dEdxView = 0;

  for (Int_t ip = 0; ip < 6; ip++) {
    // X view
    dEdxView = trkTrack->GetDEDX(ip, 0);
    if (dEdxView > 0 && trkTrack->XGood(ip)) {
      totReleaseX += dEdxView;
      nX++;
    }
    // Y view
    dEdxView = trkTrack->GetDEDX(ip, 1);
    if (dEdxView > 0 && trkTrack->YGood(ip)) {
      totReleaseY += dEdxView;
      nY++;
    }
  }

  Float_t dEdxX = totReleaseX / nX; // nX assumed > 0
  Float_t dEdxY = totReleaseY / nY; // nY assumed > 0
  
  Float_t denHigh = pow(rigMod, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
  if (dEdxX > 3.7 + 4.6 / denHigh) {
    return 0;
  }
  if (dEdxY > 3.3 + 4.9 / denHigh) {
    return 0;
  }
   
  return CUTOK;

}
1	pam-fi	1.1	/*
2			* TrkDedxHCut.cpp
3			*
4	pam-fi	1.2	* Created on: 15-mar-2010
5			* Author: S. Ricciarini
6	pam-fi	1.1	*/
7
8			/! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file /
9
10			#include "TrkDedxHCut.h"
11
12			const Double_t TrkDedxHCut::_xLow[] = { 0.126414, 0.140142, 0.152714, 0.165854, 0.178357, 0.194244, 0.211669, 0.242859,
13			0.264646, 0.293385, 0.314258, 0.360565, 0.392911, 0.450807, 0.508423, 0.603735, 0.680897, 0.754838, 0.851312,
14			1.00032, 1.36201, 1.37, 1000. }; // _nLow=23 elements
15
16			const Double_t TrkDedxHCut::_yLow[] = { 18.9221, 16.0752, 14.2214, 12.6324, 11.2421, 9.71928, 8.1965, 6.93856, 6.14407,
17			5.15095, 4.6875, 3.8268, 3.49576, 2.89989, 2.56886, 2.17161, 1.90678, 1.70816, 1.50953, 1.27754, 0.715042, 0, 0 }; // _nLow=23 elements
18
19	pam-fi	1.2	Int_t TrkDedxHCut::Check(PamLevel2 *event) {
20	pam-fi	1.1
21	pam-fi	1.3	// TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
22			if(event->GetNTracks(_trkAlg)==0)return 0;
23			ExtTrack *trkTrack = event->GetTrack(0,_trkAlg)->GetExtTrack();
24	pam-fi	1.1
25			Float_t rigMod = trkTrack->GetRigidity();
26			Float_t dedx = trkTrack->GetDEDX();
27
28	pam-fi	1.2	if (rigMod < 0.126414)
29	pam-fi	1.1	return 0; // a proton below 126.414 MV is rejected
30
31	pam-fi	1.2	// lower limit taken from N. De Simone (interpolation of _xLow)
32
33			for (Int_t i = 0; i < _nLow - 1; i++) {
34	pam-fi	1.1	if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
35	pam-fi	1.2	Double_t slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
36			Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);
37	pam-fi	1.1
38			if (dedx < y)
39			return 0;
40			}
41			}
42
43	pam-fi	1.2	// higher limit taken from N. Mori lower limit for Helium
44	pam-fi	1.1
45	pam-fi	1.2	// Compute dE/dx for X and Y separately
46			Float_t totReleaseX = 0.;
47			Float_t totReleaseY = 0.;
48			Int_t nX = 0, nY = 0;
49			Float_t dEdxView = 0;
50
51			for (Int_t ip = 0; ip < 6; ip++) {
52			// X view
53			dEdxView = trkTrack->GetDEDX(ip, 0);
54			if (dEdxView > 0 && trkTrack->XGood(ip)) {
55			totReleaseX += dEdxView;
56			nX++;
57			}
58			// Y view
59			dEdxView = trkTrack->GetDEDX(ip, 1);
60			if (dEdxView > 0 && trkTrack->YGood(ip)) {
61			totReleaseY += dEdxView;
62			nY++;
63	pam-fi	1.1	}
64			}
65
66	pam-fi	1.2	Float_t dEdxX = totReleaseX / nX; // nX assumed > 0
67			Float_t dEdxY = totReleaseY / nY; // nY assumed > 0
68
69			Float_t denHigh = pow(rigMod, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
70			if (dEdxX > 3.7 + 4.6 / denHigh) {
71			return 0;
72			}
73			if (dEdxY > 3.3 + 4.9 / denHigh) {
74			return 0;
75	pam-fi	1.1	}
76	pam-fi	1.2
77	pam-fi	1.1	return CUTOK;
78
79			}