/[PAMELA software]/PamCut/TrkCuts/TrkDedxHCut/TrkDedxHCut.cpp
ViewVC logotype

Diff of /PamCut/TrkCuts/TrkDedxHCut/TrkDedxHCut.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.1 by pam-fi, Wed May 27 13:30:08 2009 UTC revision 1.2 by pam-fi, Thu Mar 18 14:46:11 2010 UTC
# Line 1  Line 1 
1  /*  /*
2   * TrkDedxHCut.cpp   * TrkDedxHCut.cpp
3   *   *
4   *  Created on: 26-mar-2009   *  Created on: 15-mar-2010
5   *      Author: Nico De Simone, S. Ricciarini   *      Author: S. Ricciarini
6   */   */
7    
8  /*! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file */  /*! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file */
# Line 16  const Double_t TrkDedxHCut::_xLow[] = { Line 16  const Double_t TrkDedxHCut::_xLow[] = {
16  const Double_t TrkDedxHCut::_yLow[] = { 18.9221, 16.0752, 14.2214, 12.6324, 11.2421, 9.71928, 8.1965, 6.93856, 6.14407,  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      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  const Double_t TrkDedxHCut::_xHigh[] = { 0.126414, 0.392911, 0.474654, 0.508423, 0.563635, 0.680897, 0.76792, 0.851312,  Int_t TrkDedxHCut::Check(PamLevel2 *event) {
     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) {  
20    
21    TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();    TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
22    
23    Float_t rigMod = trkTrack->GetRigidity();    Float_t rigMod = trkTrack->GetRigidity();
24    Float_t dedx = trkTrack->GetDEDX();    Float_t dedx = trkTrack->GetDEDX();
25    
26  #ifdef DEBUGPAMCUT    if (rigMod < 0.126414)
   
   h_trk_he_dedx_rigmod[0]->Fill(rigMod, dedx);  
   h_trk_le_dedx_rigmod[0]->Fill(rigMod, dedx);  
   
 #endif  
   
   if (rigMod < _xHigh[0])  
27      return 0; // a proton below 126.414 MV is rejected      return 0; // a proton below 126.414 MV is rejected
28    
29    // interpolation of _xLow    // lower limit taken from N. De Simone (interpolation of _xLow)
30    for (int i = 0; i < _nLow - 1; i++) {    
31      for (Int_t i = 0; i < _nLow - 1; i++) {
32      if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {      if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
33        double slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);        Double_t slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
34        double y = _yLow[i] + slope * (rigMod - _xLow[i]);        Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);
35    
36        if (dedx < y)        if (dedx < y)
37          return 0;          return 0;
38      }      }
39    }    }
40    
41    // interpolation of _xHigh    // higher limit taken from N. Mori lower limit for Helium
   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]);  
42    
43        if (dedx > y)    // Compute dE/dx for X and Y separately
44          return 0;    Float_t totReleaseX = 0.;
45      Float_t totReleaseY = 0.;
46      Int_t nX = 0, nY = 0;
47      Float_t dEdxView = 0;
48    
49      for (Int_t ip = 0; ip < 6; ip++) {
50        // X view
51        dEdxView = trkTrack->GetDEDX(ip, 0);
52        if (dEdxView > 0 && trkTrack->XGood(ip)) {
53          totReleaseX += dEdxView;
54          nX++;
55        }
56        // Y view
57        dEdxView = trkTrack->GetDEDX(ip, 1);
58        if (dEdxView > 0 && trkTrack->YGood(ip)) {
59          totReleaseY += dEdxView;
60          nY++;
61      }      }
62    }    }
63    
64    // extrapolation above 1 TV    Float_t dEdxX = totReleaseX / nX; // nX assumed > 0
65    if (rigMod >= _xHigh[_nHigh - 1]) {    Float_t dEdxY = totReleaseY / nY; // nY assumed > 0
66      double slope = (_yHigh[_nHigh - 1] - _yHigh[_nHigh - 2]) / (_xHigh[_nHigh - 1] - _xHigh[_nHigh - 2]);    
67      double y = _yHigh[_nHigh - 1] + slope * (rigMod - _xHigh[_nHigh - 1]);    Float_t denHigh = pow(rigMod, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
68      if (dEdxX > 3.7 + 4.6 / denHigh) {
69      if (dedx > y)      return 0;
       return 0;  
70    }    }
71      if (dEdxY > 3.3 + 4.9 / denHigh) {
72  #ifdef DEBUGPAMCUT      return 0;
73      }
74    h_trk_he_dedx_rigmod[1]->Fill(rigMod, dedx);    
   h_trk_le_dedx_rigmod[1]->Fill(rigMod, dedx);  
   
 #endif  
   
75    return CUTOK;    return CUTOK;
76    
77  }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23