source: trunk/MagicSoft/Mars/mfilter/MFGeomag.cc@ 2510

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): R.K.Bock 11/2003     <mailto:rkb@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//   MFGeomag
//
//  A filter to reject Monte Carlo events based on phi/theta/charge of the
//  incident particle. Uses tables calculated by Adrian Biland, which contain
//  three parameters, used with rigidity (= particle momentum / charge) :
//         rig < min_rig:   reject unconditionally
//         rig > max_rig:   accept unconditionally
//         rig in between:  reject it with 'probability'
//  the two tables, for + and - rigidity, are stored in ASCII form in mfilter/
//
/////////////////////////////////////////////////////////////////////////////

#include "MFGeomag.h"

#include "fstream"        //for ifstream
#include "TRandom.h"      //for gRandom

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

#include "MMcEvt.hxx"


ClassImp(MFGeomag);

using namespace std;

// --------------------------------------------------------------------------
//
MFGeomag::MFGeomag(const char *cname, const char type, const Int_t val,
                           const char *name, const char *title) : fMcEvt(NULL)
{
    fContName = cname;
    Init(type, val, name, title);
}

// --------------------------------------------------------------------------
//
void MFGeomag::Init(const char type, const Int_t val,
                        const char *name, const char *title)

{
    fName  = name  ? name  : "MFGeomag";
    fTitle = title ? title : "Filter using geomagnetic field";

    fGamm_elec = kFALSE;  // logical variable, will not take gammas as electrons (default)

    AddToBranchList(Form("%s.fPartId", (const char*)fContName));
}
// --------------------------------------------------------------------------
//
Int_t MFGeomag::PreProcess(MParList *pList)
{
    //  reading of tables (variables defined as 'private')

    Float_t azim  [2*1152];      // (these variables not used)
    Float_t thet  [2*1152];
    ifstream geomagp;
    geomagp.open("mfilter/gcplus.txt");
    if (!geomagp) {
        *fLog << err <<" ERROR gcplus.txt file not found by Geomag"<<endl;
        return kFALSE;
    }
    for (int i=0; i<1152; i++)
    {
        geomagp >>azim[i]>>thet[i];
        geomagp >>fRigMin[i]>>fRigMax[i]>>fProb[i];
    }
    *fLog << inf << endl;
    *fLog << "gcplus.txt  read, first line: ";
    *fLog << Form ("FRigMin=%8f  fRigMax=%8f  fProb=%8f",
                   fRigMin[0], fRigMax[0], fProb[0]) << endl;
    geomagp.close();
    geomagp.clear();
    ifstream geomagm;
    geomagm.open("mfilter/gcminus.txt");
    if (!geomagm) {
        *fLog << err <<" ERROR gcminus.txt file not found by Geomag"<<endl;
        return kFALSE;
    }
    for (int i=0; i<1152; i++)
    {
        geomagm >>azim[i+1152]>>thet[i+1152];
        geomagm >>fRigMin[i+1152]>>fRigMax[i+1152]>>fProb[i+1152];
    }
    *fLog << "gcminus.txt read, first line: ";
    *fLog << Form ("fRigMin=%8f  fRigMax=%8f  fProb=%8f",
                   fRigMin[1152], fRigMax[1152], fProb[1152]) << endl;
    geomagm.close();
    geomagm.clear();
    //
    if (fMcEvt)
        return kTRUE;

    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    if (fMcEvt)
        return kTRUE;

    *fLog << err << dbginf << "  [MMcEvt] not found... aborting." << endl;
    return kFALSE;
}
// --------------------------------------------------------------------------
//
void MFGeomag::SetGammElec()
{
    fGamm_elec = kTRUE;  // logical variable, will take gammas as electrons
    *fLog <<" MFGeomag called to treat gammas as electrons" << endl;
    return;
}

// --------------------------------------------------------------------------
//
Int_t MFGeomag::Process()
   {
    const Float_t en =  fMcEvt->GetEnergy();       // for rigidity (set P = E)
    float 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)
    switch (fMcEvt->GetPartId())
    {
    case kGAMMA:
        if (!fGamm_elec)         //accept gammas if not set to electrons
            return kTRUE;
        indadd = 1152;           //second part of table (negative particles)
        break;

    case kHELIUM:
        rig /= 2;                //double charge
        break;

    case kPROTON:                //protons
        break;

    default:
        Int_t id = fMcEvt->GetPartId();
        *fLog << err <<" Unknown Monte Carlo Particle Id#: "<< id << endl;
        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=1;        // reject
          return kTRUE;
      }
      float r2=fRigMax[72*it+ia+indadd];
      if (rig>=r2) {
          fResult=0;        // accept
          return kTRUE;
      }
      float R = gRandom->Rndm(0);        //accept if above intermediate threshold
      float pr=fProb  [72*it+ia+indadd];    
      fResult = 0;
      if (rig < 0.5/pr*R*(r2-r1) + r1)  fResult = 1; // pretty good approximation
      return kTRUE;
   }