--- PamCut/TrkCuts/TrkDedxHCut/TrkDedxHCut.cpp	2009/05/27 13:30:08	1.1
+++ PamCut/TrkCuts/TrkDedxHCut/TrkDedxHCut.cpp	2010/03/18 14:46:11	1.2
@@ -1,8 +1,8 @@
 /*
  * TrkDedxHCut.cpp
  *
- *  Created on: 26-mar-2009
- *      Author: Nico De Simone, S. Ricciarini
+ *  Created on: 15-mar-2010
+ *      Author: S. Ricciarini
  */
 
 /*! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file */
@@ -16,68 +16,62 @@
 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
 
-const Double_t TrkDedxHCut::_xHigh[] = { 0.126414, 0.392911, 0.474654, 0.508423, 0.563635, 0.680897, 0.76792, 0.851312,
-    0.960115, 1.12068, 1.22122, 1.50085, 1.69267, 1.87649, 2.1903, 2.65934, 3.75055, 6.21589, 11.0139, 24.2773, 1000. }; // _nHigh=21 elements
-
-const Double_t TrkDedxHCut::_yHigh[] = { 24.3512, 22.696, 21.107, 20.0477, 18.5249, 15.8766, 13.6917, 11.5731, 10.1827,
-    8.26271, 7.2696, 5.5482, 5.01854, 4.6875, 4.15784, 3.83697, 3.59789, 3.42119, 3.36335, 3.36335, 3.36335 }; // _nHigh=21 elements
-
-
-int TrkDedxHCut::Check(PamLevel2 *event) {
+Int_t TrkDedxHCut::Check(PamLevel2 *event) {
 
   TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
 
   Float_t rigMod = trkTrack->GetRigidity();
   Float_t dedx = trkTrack->GetDEDX();
 
-#ifdef DEBUGPAMCUT
-
-  h_trk_he_dedx_rigmod[0]->Fill(rigMod, dedx);
-  h_trk_le_dedx_rigmod[0]->Fill(rigMod, dedx);
-
-#endif
-
-  if (rigMod < _xHigh[0])
+  if (rigMod < 0.126414)
     return 0; // a proton below 126.414 MV is rejected
 
-  // interpolation of _xLow
-  for (int i = 0; i < _nLow - 1; i++) {
+  // 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 slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
-      double y = _yLow[i] + slope * (rigMod - _xLow[i]);
+      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;
     }
   }
 
-  // interpolation of _xHigh
-  for (int i = 0; i < _nHigh - 1; i++) {
-    if (_xHigh[i] <= rigMod && rigMod < _xHigh[i + 1]) {
-      double slope = (_yHigh[i + 1] - _yHigh[i]) / (_xHigh[i + 1] - _xHigh[i]);
-      double y = _yHigh[i] + slope * (rigMod - _xHigh[i]);
+  // higher limit taken from N. Mori lower limit for Helium
 
-      if (dedx > y)
-        return 0;
+  // 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++;
     }
   }
 
-  // extrapolation above 1 TV
-  if (rigMod >= _xHigh[_nHigh - 1]) {
-    double slope = (_yHigh[_nHigh - 1] - _yHigh[_nHigh - 2]) / (_xHigh[_nHigh - 1] - _xHigh[_nHigh - 2]);
-    double y = _yHigh[_nHigh - 1] + slope * (rigMod - _xHigh[_nHigh - 1]);
-
-    if (dedx > y)
-      return 0;
+  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;
   }
-
-#ifdef DEBUGPAMCUT
-
-  h_trk_he_dedx_rigmod[1]->Fill(rigMod, dedx);
-  h_trk_le_dedx_rigmod[1]->Fill(rigMod, dedx);
-
-#endif
-
+  if (dEdxY > 3.3 + 4.9 / denHigh) {
+    return 0;
+  }
+   
   return CUTOK;
 
 }