--- PamCut/TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp	2009/06/10 12:37:57	1.1
+++ PamCut/TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp	2010/05/03 14:36:27	1.2
@@ -1,8 +1,8 @@
 /*
  * TrkDedxHeCut.cpp
  *
- *  Created on: 19-may-2009
- *      Author: N. Nikonov
+ *  Created on: 28/gen/2010
+ *      Author: Nicola Mori
  */
 
 /*! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. */
@@ -11,57 +11,91 @@
 
 int TrkDedxHeCut::Check(PamLevel2 *event) {
 
-  TrkTrack* trk2 = event->GetTrack(0)->GetTrkTrack();
-  Double_t rig = trk2->GetRigidity();
-  Double_t dedxtr = trk2->GetDEDX();
-
-#ifdef DEBUGPAMCUT
-  h_dedx_rig_before->Fill(rig, dedxtr);
-#endif
-
-  if ( !CutHeDEDXrig(rig,dedxtr,1.15,_cuthededxrig) ||
-       CutHeDEDXrig(rig,dedxtr,0.5,_cuthededxrig)) return 0;
-
-#ifdef DEBUGPAMCUT
-  h_dedx_rig_after->Fill(rig, dedxtr);
-#endif
-
-  return CUTOK;
-}
-
-
-TF1* TrkDedxHeCut::BBbetagamma(Float_t m_GeV, Float_t z, Int_t Z, Float_t A, Float_t I_eV, Double_t xmin, Double_t xmax, Float_t C_positive, Float_t a, Float_t m, Float_t x0, Float_t x1, Float_t scale)
-{
-  TF1* f;
-  if(C_positive != 0.) {
-    // do density correction
-    // x<x_0 ( delta = 0)
-    // I don't need a function in this range
-    // x0<x<x1
-    // I build the function as a string.
-    TString delta2 = "(4.6052*log([1]*x/[0])/log(10)-[5]+[6]*pow(([9]-log([1]*x/[0])/log(10)),[7])) * (([8] <= log([1]*x/[0])) && ( log([1]*x/[0])/log(10)  < [9] ) )";
-    // x>x1
-    TString delta3 = "(4.6052*log([1]*x/[0])/log(10)-[5]) * (log([1]*x/[0])/log(10) >= [9])";
-
-    TString fstring = "[10]*(.3071 * ([2]/[3]) * [1]^2 / ([0]^2/([1]*x)^2+1)^-1 * ( 1./2 * log( (2 * 511e3 * ([1]*x/[0])^2 )^2 ) - log([4]) - ([0]^2/([1]*x)^2+1)^-1 - (";
-    fstring += delta2.Data();
-    fstring += "+";
-    fstring += delta3.Data();
-    fstring += ")/2))";
-    f = new TF1("bethebloch", fstring.Data());
+  TrkTrack *track = event->GetTrack(0)->GetTrkTrack();
+  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;
+  /*
+   ToFLevel2 *tofL2 = event->GetToFLevel2();
+   int trkSeqNo = 0;
+   for (int i = 0; i < event->GetToFLevel2()->ntrk(); i++) {
+   if (tofL2->GetToFTrkVar(i)->trkseqno == track->GetSeqNo()) {
+   trkSeqNo = i;
+   break;
+   }
+   }
+   if ((tofL2->GetdEdx(trkSeqNo, 0, 100) + tofL2->GetdEdx(trkSeqNo, 1, 100) + tofL2->GetdEdx(trkSeqNo, 2, 100)
+   + tofL2->GetdEdx(trkSeqNo, 3, 100)) / 4. < 2) {
+   cout << "Check" << endl;
+   for (int ip = 0; ip < 6; ip++) {
+   cout << "TOF" << ip << ": " << tofL2->GetdEdx(trkSeqNo, ip, 100) << endl;
+   dEdxView = track->GetDEDX(ip, 0);
+   if (dEdxView > 0 && track->XGood(ip))
+   cout << "X" << ip << ": " << dEdxView << endl;
+   dEdxView = track->GetDEDX(ip, 1);
+   if (dEdxView > 0 && track->YGood(ip))
+   cout << "Y" << ip << ": " << dEdxView << endl;
+   }
+   }
+   */
+  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;
+    }
   }
-  else
-    f = new TF1("bethebloch", ".3071 * ([2]/[3]) * [1]^2 / ([0]^2/([1]*x)^2+1)^-1 * ( 1./2 * log( (2 * 511e3 * ([1]*x/[0])^2 )^2 ) - log([4]) - ([0]^2/([1]*x)^2+1)^-1)");
-  f->SetParameters(m_GeV, z, Z, A, I_eV, C_positive, a, m, x0, x1,scale);
 
-  return f;
-}
+  // 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;*/
 
-Bool_t TrkDedxHeCut::CutHeDEDXrig(Float_t rig, Float_t dedx,
-					 Float_t factor, TF1* _cuthededxrig){
-  if(factor != 1.) dedx = dedx*factor;
-  if(rig < 0.2)  return kFALSE;
-  if(dedx > 0.55*_cuthededxrig->Eval(rig)) return kTRUE;
-  else return kFALSE;
+  return CUTOK;
 }