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	/*
2	* TrkDedxHeCut.cpp
3	*
4	* Created on: 28/gen/2010
5	* Author: Nicola Mori
6	*/
7
8	/! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. /
9
10	#include "TrkDedxHeCut.h"
11
12	int TrkDedxHeCut::Check(PamLevel2 *event) {
13
14	// 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	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
47	// 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	}
56
57	// 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
81	return CUTOK;
82	}