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

Diff of /PamCut/TrkCuts/TrkDedxHeCut/TrkDedxHeCut.cpp

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

revision 1.1 by pam-fi, Wed Jun 10 12:37:57 2009 UTC revision 1.3 by pam-fi, Wed Mar 25 17:38:08 2015 UTC
# Line 1  Line 1 
1  /*  /*
2   * TrkDedxHeCut.cpp   * TrkDedxHeCut.cpp
3   *   *
4   *  Created on: 19-may-2009   *  Created on: 28/gen/2010
5   *      Author: N. Nikonov   *      Author: Nicola Mori
6   */   */
7    
8  /*! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. */  /*! @file TrkDedxHeCut.cpp The TrkDedxHeCut class implementation file. */
# Line 11  Line 11 
11    
12  int TrkDedxHeCut::Check(PamLevel2 *event) {  int TrkDedxHeCut::Check(PamLevel2 *event) {
13    
14    TrkTrack* trk2 = event->GetTrack(0)->GetTrkTrack();  //  TrkTrack *track = event->GetTrack(0)->GetTrkTrack();
15    Double_t rig = trk2->GetRigidity();      if(event->GetNTracks(_trkAlg)==0)return 0;    
16    Double_t dedxtr = trk2->GetDEDX();      ExtTrack *track = event->GetTrack(0,_trkAlg)->GetExtTrack();
17        
18  #ifdef DEBUGPAMCUT  
19    h_dedx_rig_before->Fill(rig, dedxtr);    float R = track->GetRigidity();
20  #endif    //float dEdx = track->GetDEDX();
21    
22    if ( !CutHeDEDXrig(rig,dedxtr,1.15,_cuthededxrig) ||    // New smart dE/dx computation
23         CutHeDEDXrig(rig,dedxtr,0.5,_cuthededxrig)) return 0;    // This procedure removes high releases and computes dE/dx
24      float totReleaseX = 0., maxReleaseX = 0.;
25  #ifdef DEBUGPAMCUT    float totReleaseY = 0., maxReleaseY = 0.;
26    h_dedx_rig_after->Fill(rig, dedxtr);    int nX = 0, nY = 0;
27  #endif    float dEdxView = 0;
28      for (int ip = 0; ip < 6; ip++) {
29    return CUTOK;      // X view
30  }      dEdxView = track->GetDEDX(ip, 0);
31        if (dEdxView > 0 && track->XGood(ip)) {
32          totReleaseX += dEdxView;
33  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)        nX++;
34  {        if (dEdxView > maxReleaseX)
35    TF1* f;          maxReleaseX = dEdxView;
36    if(C_positive != 0.) {      }
37      // do density correction      // Y view
38      // x<x_0 ( delta = 0)      dEdxView = track->GetDEDX(ip, 1);
39      // I don't need a function in this range      if (dEdxView > 0 && track->YGood(ip)) {
40      // x0<x<x1        totReleaseY += dEdxView;
41      // I build the function as a string.        nY++;
42      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] ) )";        if (dEdxView > maxReleaseY)
43      // x>x1          maxReleaseY = dEdxView;
44      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());  
45    }    }
   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);  
46    
47    return f;    // Discard highest release, eventually
48  }    if (maxReleaseX > -9. + 4. * totReleaseX / (float) nX) {
49        totReleaseX -= maxReleaseX;
50        nX--;
51      }
52      if (maxReleaseY > -8. + 4. * totReleaseY / (float) nY) {
53        totReleaseY -= maxReleaseY;
54        nY--;
55      }
56    
57      // Compute dE/dx
58      //float dEdx = (totReleaseX + totReleaseY) / (float) (nX + nY);
59      float dEdxX = totReleaseX / nX;
60      float dEdxY = totReleaseY / nY;
61    
62      // Analyze the event
63      // X
64      float denLow = pow(R, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
65      float denHigh = pow(R, 1.5);
66      if (dEdxX < 3.7 + 4.6 / denLow)
67        return 0;
68      if (dEdxX > 8.9 + 17. / denHigh)
69        return 0;
70      // Y
71      if (dEdxY < 3.3 + 4.9 / denLow)
72        return 0;
73      if (dEdxY > 8.0 + 17. / denHigh)
74        return 0;
75      //X+Y
76      /*if (dEdx < 3.7 + 4.8 / (R * R))
77        return 0;
78      if (dEdx > 8.9 + 17. / (R * R))
79        return 0;*/
80    
81  Bool_t TrkDedxHeCut::CutHeDEDXrig(Float_t rig, Float_t dedx,    return CUTOK;
                                          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;  
82  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23