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

Diff of /PamCut/TrkCuts/TrkChi2DeflTimeCut/TrkChi2DeflTimeCut.cpp

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

revision 1.2 by pam-fi, Thu Nov 12 16:41:43 2009 UTC revision 1.6 by pam-fi, Fri Sep 2 16:28:17 2011 UTC
# Line 7  Line 7 
7    
8  #include "TrkChi2DeflTimeCut.h"  #include "TrkChi2DeflTimeCut.h"
9    
10  double TrkChi2DeflTimeCut::_GetChi2(PamLevel2 *event){  double TrkChi2DeflTimeCut::_GetChi2(PamLevel2 *event) {
11    return event->GetTrack(0)->GetTrkTrack()->chi2;    return event->GetTrack(0)->GetTrkTrack()->chi2;
12  }  }
13    
# Line 16  int TrkChi2DeflTimeCut::Check(PamLevel2 Line 16  int TrkChi2DeflTimeCut::Check(PamLevel2
16    double chi2 = _GetChi2(event);    double chi2 = _GetChi2(event);
17    if (chi2 < 0.)    if (chi2 < 0.)
18      return 0;      return 0;
19    Double_t eta2 = event->GetTrack(0)->GetTrkTrack()->GetDeflection(); //Sign doesn't matter    Double_t etaMod = fabs(event->GetTrack(0)->GetTrkTrack()->GetDeflection());
20    eta2 *= eta2;  // Now eta2 is deflection^2    Int_t nHitX = event->GetTrack(0)->GetTrkTrack()->GetNX();
21      
22  #ifdef DEBUGPAMCUT  //  cout << " nHitX " << nHitX << endl;
23      
24    h_trk_chi2_defl[0]->Fill(event->GetTrack(0)->GetTrkTrack()->GetDeflection(), pow(chi2,0.25));    if (nHitX == 3) {
25        for (UInt_t i = 0; i < _iDayFirst_nHitX3.size(); i++) {
26  #endif        _iDayFirst.push_back(_iDayFirst_nHitX3[i]);
27          _iDayLast.push_back(_iDayLast_nHitX3[i]);
28          for (int iPar=0; iPar<_nPar; iPar++) {
29            _p[iPar].push_back(_p_nHitX3[iPar][i]);
30          }
31        }  
32      }
33      else if (nHitX >= 4) {
34        for (UInt_t i = 0; i < _iDayFirst_nHitX4.size(); i++) {
35          _iDayFirst.push_back(_iDayFirst_nHitX4[i]);
36          _iDayLast.push_back(_iDayLast_nHitX4[i]);
37          for (int iPar=0; iPar<_nPar; iPar++) {
38            _p[iPar].push_back(_p_nHitX4[iPar][i]);
39          }
40        }  
41      }
42      else {
43        return 0;
44      }
45    
46    _time.Set(event->GetOrbitalInfo()->absTime, kFALSE, 0, kFALSE);    _time.Set(event->GetOrbitalInfo()->absTime, kFALSE, 0, kFALSE);
47    // TTimestamp::GetDate() returns date in format YYYYMMDD so to compare it    // TTimestamp::GetDate() returns date in format YYYYMMDD so to compare it
48    // with 4-digits dates (YYMM) we must subtract 20000000 and then divide by 100 (divide by integer)    // with 6-digits dates (YYMMDD) we must subtract 20000000
49    
50    Int_t iMonthSel = (_time.GetDate(kFALSE) - 20000000)/100; // note: GetDate returns an unsigned integer    Int_t iDaySel = (_time.GetDate(kFALSE) - 20000000); // note: GetDate returns an unsigned integer
51        //  cout << "day " << iDaySel << endl;
52    for (UInt_t i=0; i<_iMonth.size(); i++) {  
53      if (_iMonth[i]==iMonthSel) {    bool found = false;
54        _p0sel=_p0[i];    for (UInt_t i = 0; i < _iDayFirst.size(); i++) {
55        _p1sel=_p1[i];      if (_iDayFirst[i] <= iDaySel && iDaySel <= _iDayLast[i]) {
56        _p2sel=_p2[i];        for (int iPar=0; iPar<_nPar; iPar++) {
57  //      cout << "FOUND: month " << iMonthSel << endl;          _pSel[iPar] = _p[iPar][i];
58  //      cout << _p0sel << " " << _p1sel << " " << _p2sel << endl;        }
59          
60    //      cout << "FOUND: day " << iDaySel << endl;
61    //      cout << _iDayFirst[i] << " " << _iDayLast[i] << endl;
62    //      cout << _pSel[0] << " " << _pSel[1] << " " << _pSel[2] << " " << _pSel[3] << endl;
63    
64          found = true;
65        break;        break;
66      }      }
67    }    }
68      
69    if (chi2 > _p0sel + _p1sel * eta2 + _p2sel * eta2 * eta2) {    if (!found) {
70        cout << "not found" << endl;
71        cout << endl;
72      return 0;      return 0;
73    }    }
74    
75  #ifdef DEBUGPAMCUT  //  for (int iPar=0; iPar<_nPar; iPar++) {
76    //    cout << " p" << iPar << ": " << _pSel[iPar] << endl;
77    //  }
78    //  cout << endl;
79    
80    h_trk_chi2_defl[1]->Fill(event->GetTrack(0)->GetTrkTrack()->GetDeflection(), pow(chi2,0.25));    double chi2max;
81      
82      if (_nPar==4) {
83      
84        chi2max=_pSel[0] + _pSel[1] * pow(etaMod, _pSel[2]) * (1. + pow(_pSel[3] * etaMod, 2));
85    
86    //    cout << "chi2max " << chi2max << endl;
87    
88      }
89      else if (_nPar==5) {
90        chi2max=_pSel[0] + _pSel[1] * pow(etaMod, _pSel[2]) * (1. + pow(_pSel[3] * etaMod, _pSel[4]));
91      }
92      else {
93        chi2max=-1;
94      }
95    
96    //  cout << "chi2 " << chi2 << endl;
97      
98      if (chi2 > chi2max) {
99      
100    //    cout << "KO" << endl;
101        
102        return 0;
103    
104      }
105    
106  #endif  //  cout << "OK" << endl;
107    
108    return CUTOK;    return CUTOK;
109    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23