Changeset 2990


Ignore:
Timestamp:
01/30/04 15:28:12 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r2984 r2990  
    4646   * mraw/MRawRunHeader.[h,cc]:
    4747     - initialize fRunType to new enum kRTNone
     48
     49   * mfilter/MFGeomag.cc:
     50     - set fResult to kFALSE at the beginning of Process
     51     - according to this removed setting fResult=kFALSE before return
     52     - replaced some float by Float_t
     53     - added some const-qualifiers
     54     - replaced (rig<0.5/pr*R*(r2-r1)+r1) by (rig-r1)*pr<rnd
     55       with rnd = R * (r2-r1)/2 to make sure that we cannot
     56       devide by 0
    4857
    4958
  • trunk/MagicSoft/Mars/mfilter/MFGeomag.cc

    r2802 r2990  
    145145//
    146146Int_t MFGeomag::Process()
    147    {
     147{
     148    fResult = kFALSE;
     149
    148150    const Float_t en =  fMcEvt->GetEnergy();       // for rigidity (set P = E)
    149     float rig = en;
     151    Float_t rig = en;
    150152    const Float_t az =  fMcEvt->GetTelescopePhi(); // charge theta phi are entries in table
    151153    const Float_t th =  fMcEvt->GetTelescopeTheta();
    152154
    153     Int_t indadd=0;         //first part of table (positive particles)
     155    Int_t indadd=0;              //first part of table (positive particles)
    154156    switch (fMcEvt->GetPartId())
    155157    {
    156158    case kGAMMA:
    157         if (!fGammaElectron)         //accept gammas if not set to electrons
    158           {
    159             fResult = kFALSE;
    160             return kTRUE;
    161           }
     159        if (!fGammaElectron)     //accept gammas if not set to electrons
     160            return kTRUE;
    162161        indadd = 1152;           //second part of table (negative particles)
    163162        break;
     
    168167
    169168    case kPROTON:                //protons
     169    case kPOSITRON:              //positrons
    170170        break;
    171171
    172172    case kELECTRON:              //electrons
    173173        indadd = 1152;           //second part of table (negative particles)
    174         break;         
    175 
    176     case kPOSITRON:              //positrons
    177174        break;
    178175
     
    181178        return kFALSE;
    182179    }
    183              // here is the cut for charged particles using the table
    184       int it=(int)(th*11.459156);    // table steps are in 5 deg = 1/11.459 rads
    185       int ia=(int)(az*11.459156);
    186       ia = (ia+36) % 72;             // azimuth definitions differ by 180 deg
    187 
    188       float r1=fRigMin[72*it+ia+indadd];
    189       if (rig<=r1) {
    190           fResult=kTRUE;    // reject
    191           return kTRUE;
    192       }
    193 
    194       float r2=fRigMax[72*it+ia+indadd];
    195       if (rig>=r2) {
    196           fResult=kFALSE;   // accept
    197           return kTRUE;
    198       }
    199 
    200       float R = gRandom->Rndm(0);        // accept if above intermediate threshold
    201       float pr=fProb[72*it+ia+indadd];
    202       fResult = kFALSE;
    203 
    204       if (rig < 0.5/pr*R*(r2-r1) + r1)
    205           fResult = kTRUE;               // pretty good approximation
    206 
    207       return kTRUE;
    208    }
     180
     181    // here is the cut for charged particles using the table
     182
     183    int it=(int)(th*11.459156);    // table steps are in 5 deg = 1/11.459 rads
     184    int ia=(int)(az*11.459156);
     185
     186    ia = (ia+36) % 72;             // azimuth definitions differ by 180 deg
     187
     188    const Float_t r1=fRigMin[72*it+ia+indadd];
     189    if (rig<=r1)
     190    {
     191        fResult=kTRUE;   // reject
     192        return kTRUE;
     193    }
     194
     195    const Float_t r2=fRigMax[72*it+ia+indadd];
     196    if (rig>=r2)
     197        return kTRUE;   // accept
     198
     199    const Float_t pr=fProb[72*it+ia+indadd];
     200
     201    // accept if above intermediate threshold
     202    const Float_t rnd = (r2-r1)/2 * gRandom->Rndm(0);
     203
     204    if ((rig-r1)*pr < rnd)
     205        fResult = kTRUE;                // pretty good approximation
     206
     207    return kTRUE;
     208}
Note: See TracChangeset for help on using the changeset viewer.