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	/*
2	* TrkDedxHCut.cpp
3	*
4	* Created on: 15-mar-2010
5	* Author: S. Ricciarini
6	*/
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	Int_t TrkDedxHCut::Check(PamLevel2 *event) {
20
21	// TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
22	if(event->GetNTracks(_trkAlg)==0)return 0;
23	ExtTrack *trkTrack = event->GetTrack(0,_trkAlg)->GetExtTrack();
24
25	Float_t rigMod = trkTrack->GetRigidity();
26	Float_t dedx = trkTrack->GetDEDX();
27
28	if (rigMod < 0.126414)
29	return 0; // a proton below 126.414 MV is rejected
30
31	// lower limit taken from N. De Simone (interpolation of _xLow)
32
33	for (Int_t i = 0; i < _nLow - 1; i++) {
34	if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
35	Double_t slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
36	Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);
37
38	if (dedx < y)
39	return 0;
40	}
41	}
42
43	// higher limit taken from N. Mori lower limit for Helium
44
45	// 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	}
64	}
65
66	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	}
76
77	return CUTOK;
78
79	}