Changeset 2990 for trunk/MagicSoft/Mars/mfilter/MFGeomag.cc

Timestamp:

01/30/04 15:28:12 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mfilter/MFGeomag.cc (modified) (3 diffs)

trunk/MagicSoft/Mars/mfilter/MFGeomag.cc

@@ -145 +145 @@
 //
 Int_t MFGeomag::Process()
-   {
+{
+    fResult = kFALSE;
+
     const Float_t en =  fMcEvt->GetEnergy();       // for rigidity (set P = E)
-    float rig = en;
+    Float_t rig = en;
     const Float_t az =  fMcEvt->GetTelescopePhi(); // charge theta phi are entries in table
     const Float_t th =  fMcEvt->GetTelescopeTheta();
 
-    Int_t indadd=0;         //first part of table (positive particles)
+    Int_t indadd=0;              //first part of table (positive particles)
     switch (fMcEvt->GetPartId())
     {
     case kGAMMA:
-        if (!fGammaElectron)         //accept gammas if not set to electrons
-          {
-            fResult = kFALSE;
-            return kTRUE;
-          }
+        if (!fGammaElectron)     //accept gammas if not set to electrons
+            return kTRUE;
         indadd = 1152;           //second part of table (negative particles)
         break;
@@ -168 +167 @@
 
     case kPROTON:                //protons
+    case kPOSITRON:              //positrons
         break;
 
     case kELECTRON:              //electrons
         indadd = 1152;           //second part of table (negative particles)
-        break;          
-
-    case kPOSITRON:              //positrons
         break;
 
@@ -181 +178 @@
         return kFALSE;
     }
-             // here is the cut for charged particles using the table
-      int it=(int)(th*11.459156);    // table steps are in 5 deg = 1/11.459 rads
-      int ia=(int)(az*11.459156);
-      ia = (ia+36) % 72;             // azimuth definitions differ by 180 deg
-
-      float r1=fRigMin[72*it+ia+indadd];
-      if (rig<=r1) {
-          fResult=kTRUE;    // reject
-          return kTRUE;
-      }
-
-      float r2=fRigMax[72*it+ia+indadd];
-      if (rig>=r2) {
-          fResult=kFALSE;   // accept
-          return kTRUE;
-      }
-
-      float R = gRandom->Rndm(0);        // accept if above intermediate threshold
-      float pr=fProb[72*it+ia+indadd];
-      fResult = kFALSE;
-
-      if (rig < 0.5/pr*R*(r2-r1) + r1)
-          fResult = kTRUE;               // pretty good approximation
-
-      return kTRUE;
-   }
+
+    // here is the cut for charged particles using the table
+
+    int it=(int)(th*11.459156);    // table steps are in 5 deg = 1/11.459 rads
+    int ia=(int)(az*11.459156);
+
+    ia = (ia+36) % 72;             // azimuth definitions differ by 180 deg
+
+    const Float_t r1=fRigMin[72*it+ia+indadd];
+    if (rig<=r1)
+    {
+        fResult=kTRUE;   // reject
+        return kTRUE;
+    }
+
+    const Float_t r2=fRigMax[72*it+ia+indadd];
+    if (rig>=r2)
+        return kTRUE;   // accept
+
+    const Float_t pr=fProb[72*it+ia+indadd];
+
+    // accept if above intermediate threshold
+    const Float_t rnd = (r2-r1)/2 * gRandom->Rndm(0);
+
+    if ((rig-r1)*pr < rnd)
+        fResult = kTRUE;                // pretty good approximation
+
+    return kTRUE;
+}