TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp

/*
 * TrkDedxHeCut.cpp
 *
 *  Created on: 28/gen/2010
 *      Author: Nicola Mori
 */

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

#include "TrkDedxHeCut.h"

int TrkDedxHeCut::Check(PamLevel2 *event) {

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

  float R = track->GetRigidity();
  //float dEdx = track->GetDEDX();

  // New smart dE/dx computation
  // This procedure removes high releases and computes dE/dx
  float totReleaseX = 0., maxReleaseX = 0.;
  float totReleaseY = 0., maxReleaseY = 0.;
  int nX = 0, nY = 0;
  float dEdxView = 0;
  for (int ip = 0; ip < 6; ip++) {
    // X view
    dEdxView = track->GetDEDX(ip, 0);
    if (dEdxView > 0 && track->XGood(ip)) {
      totReleaseX += dEdxView;
      nX++;
      if (dEdxView > maxReleaseX)
        maxReleaseX = dEdxView;
    }
    // Y view
    dEdxView = track->GetDEDX(ip, 1);
    if (dEdxView > 0 && track->YGood(ip)) {
      totReleaseY += dEdxView;
      nY++;
      if (dEdxView > maxReleaseY)
        maxReleaseY = dEdxView;
    }
  }

  // Discard highest release, eventually
  if (maxReleaseX > -9. + 4. * totReleaseX / (float) nX) {
    totReleaseX -= maxReleaseX;
    nX--;
  }
  if (maxReleaseY > -8. + 4. * totReleaseY / (float) nY) {
    totReleaseY -= maxReleaseY;
    nY--;
  }

  // Compute dE/dx
  //float dEdx = (totReleaseX + totReleaseY) / (float) (nX + nY);
  float dEdxX = totReleaseX / nX;
  float dEdxY = totReleaseY / nY;

  // Analyze the event
  // X
  float denLow = pow(R, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
  float denHigh = pow(R, 1.5);
  if (dEdxX < 3.7 + 4.6 / denLow)
    return 0;
  if (dEdxX > 8.9 + 17. / denHigh)
    return 0;
  // Y
  if (dEdxY < 3.3 + 4.9 / denLow)
    return 0;
  if (dEdxY > 8.0 + 17. / denHigh)
    return 0;
  //X+Y
  /*if (dEdx < 3.7 + 4.8 / (R * R))
    return 0;
  if (dEdx > 8.9 + 17. / (R * R))
    return 0;*/

  return CUTOK;
}
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.3	// TrkTrack *track = event->GetTrack(0)->GetTrkTrack();
15			if(event->GetNTracks(_trkAlg)==0)return 0;
16			ExtTrack *track = event->GetTrack(0,_trkAlg)->GetExtTrack();
17
18
19	pam-fi	1.2	float R = track->GetRigidity();
20			//float dEdx = track->GetDEDX();
21
22			// New smart dE/dx computation
23			// This procedure removes high releases and computes dE/dx
24			float totReleaseX = 0., maxReleaseX = 0.;
25			float totReleaseY = 0., maxReleaseY = 0.;
26			int nX = 0, nY = 0;
27			float dEdxView = 0;
28			for (int ip = 0; ip < 6; ip++) {
29			// X view
30			dEdxView = track->GetDEDX(ip, 0);
31			if (dEdxView > 0 && track->XGood(ip)) {
32			totReleaseX += dEdxView;
33			nX++;
34			if (dEdxView > maxReleaseX)
35			maxReleaseX = dEdxView;
36			}
37			// Y view
38			dEdxView = track->GetDEDX(ip, 1);
39			if (dEdxView > 0 && track->YGood(ip)) {
40			totReleaseY += dEdxView;
41			nY++;
42			if (dEdxView > maxReleaseY)
43			maxReleaseY = dEdxView;
44			}
45			}
46	pam-fi	1.1
47	pam-fi	1.2	// Discard highest release, eventually
48			if (maxReleaseX > -9. + 4. * totReleaseX / (float) nX) {
49			totReleaseX -= maxReleaseX;
50			nX--;
51			}
52			if (maxReleaseY > -8. + 4. * totReleaseY / (float) nY) {
53			totReleaseY -= maxReleaseY;
54			nY--;
55	pam-fi	1.1	}
56
57	pam-fi	1.2	// Compute dE/dx
58			//float dEdx = (totReleaseX + totReleaseY) / (float) (nX + nY);
59			float dEdxX = totReleaseX / nX;
60			float dEdxY = totReleaseY / nY;
61
62			// Analyze the event
63			// X
64			float denLow = pow(R, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
65			float denHigh = pow(R, 1.5);
66			if (dEdxX < 3.7 + 4.6 / denLow)
67			return 0;
68			if (dEdxX > 8.9 + 17. / denHigh)
69			return 0;
70			// Y
71			if (dEdxY < 3.3 + 4.9 / denLow)
72			return 0;
73			if (dEdxY > 8.0 + 17. / denHigh)
74			return 0;
75			//X+Y
76			/if (dEdx < 3.7 + 4.8 / (R R))
77			return 0;
78			if (dEdx > 8.9 + 17. / (R * R))
79			return 0;*/
80	pam-fi	1.1
81	pam-fi	1.2	return CUTOK;
82	pam-fi	1.1	}