/[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.3 by pam-fi, Wed Mar 25 17:38:08 2015 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  
20    
21  const Double_t TrkDedxHCut::_yHigh[] = { 24.3512, 22.696, 21.107, 20.0477, 18.5249, 15.8766, 13.6917, 11.5731, 10.1827,  //  TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
22      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      if(event->GetNTracks(_trkAlg)==0)return 0;    
23        ExtTrack *trkTrack = event->GetTrack(0,_trkAlg)->GetExtTrack();
   
 int TrkDedxHCut::Check(PamLevel2 *event) {  
   
   TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();  
24    
25    Float_t rigMod = trkTrack->GetRigidity();    Float_t rigMod = trkTrack->GetRigidity();
26    Float_t dedx = trkTrack->GetDEDX();    Float_t dedx = trkTrack->GetDEDX();
27    
28  #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])  
29      return 0; // a proton below 126.414 MV is rejected      return 0; // a proton below 126.414 MV is rejected
30    
31    // interpolation of _xLow    // lower limit taken from N. De Simone (interpolation of _xLow)
32    for (int i = 0; i < _nLow - 1; i++) {    
33      for (Int_t i = 0; i < _nLow - 1; i++) {
34      if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {      if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
35        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]);
36        double y = _yLow[i] + slope * (rigMod - _xLow[i]);        Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);
37    
38        if (dedx < y)        if (dedx < y)
39          return 0;          return 0;
40      }      }
41    }    }
42    
43    // 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]);  
44    
45        if (dedx > y)    // Compute dE/dx for X and Y separately
46          return 0;    Float_t totReleaseX = 0.;
47      Float_t totReleaseY = 0.;
48      Int_t nX = 0, nY = 0;
49      Float_t dEdxView = 0;
50    
51      for (Int_t ip = 0; ip < 6; ip++) {
52        // X view
53        dEdxView = trkTrack->GetDEDX(ip, 0);
54        if (dEdxView > 0 && trkTrack->XGood(ip)) {
55          totReleaseX += dEdxView;
56          nX++;
57        }
58        // Y view
59        dEdxView = trkTrack->GetDEDX(ip, 1);
60        if (dEdxView > 0 && trkTrack->YGood(ip)) {
61          totReleaseY += dEdxView;
62          nY++;
63      }      }
64    }    }
65    
66    // extrapolation above 1 TV    Float_t dEdxX = totReleaseX / nX; // nX assumed > 0
67    if (rigMod >= _xHigh[_nHigh - 1]) {    Float_t dEdxY = totReleaseY / nY; // nY assumed > 0
68      double slope = (_yHigh[_nHigh - 1] - _yHigh[_nHigh - 2]) / (_xHigh[_nHigh - 1] - _xHigh[_nHigh - 2]);    
69      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)
70      if (dEdxX > 3.7 + 4.6 / denHigh) {
71      if (dedx > y)      return 0;
       return 0;  
72    }    }
73      if (dEdxY > 3.3 + 4.9 / denHigh) {
74  #ifdef DEBUGPAMCUT      return 0;
75      }
76    h_trk_he_dedx_rigmod[1]->Fill(rigMod, dedx);    
   h_trk_le_dedx_rigmod[1]->Fill(rigMod, dedx);  
   
 #endif  
   
77    return CUTOK;    return CUTOK;
78    
79  }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23