KVIDINDRACsI.cpp

/***************************************************************************
                          KVIDINDRACsI.cpp  -  description
                             -------------------
    begin                : Fri Feb 20 2004
    copyright            : (C) 2004 by J.D. Frankland
    email                : frankland@ganil.fr
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "KVIDINDRACsI.h"
#include "TMath.h"
#include "KVIdentificationResult.h"
#include "TRandom.h"
#include <KVIDGCsI.h>
#include <KVINDRADetector.h>
#include <KVReconstructedNucleus.h>
 
ClassImp(KVIDINDRACsI)
 
 
 
KVIDINDRACsI::KVIDINDRACsI()
{
   CsIGrid = 0;
 
   // thresholds used by filter for Rapide-Lente identifications
   fThresMin[0][0] = 1;
   fThresMax[0][0] = 2;    // protons
   fThresMin[0][1] = 2;
   fThresMax[0][1] = 6;    // deutons
   fThresMin[0][2] = 5;
   fThresMax[0][2] = 11;   // tritons
   fThresMin[0][3] = -1;
   fThresMax[0][3] = -1;   // 4H - NO!
   fThresMin[1][0] = -1;
   fThresMax[1][0] = -1;   // 1He - NO!
   fThresMin[1][1] = -1;
   fThresMax[1][1] = -1;   // 2He - NO!
   fThresMin[1][2] = 20;
   fThresMax[1][2] = 40;   // 3He
   fThresMin[1][3] = 1;
   fThresMax[1][3] = 3;    // alphas
 
   /* in principle all CsI telescopes can identify mass & charge */
   SetHasMassID(kTRUE);
}
 
 
 
 
 
Bool_t KVIDINDRACsI::Identify(KVIdentificationResult* IDR, Double_t x, Double_t y)
{
   // Particle identification and code setting using identification grid
   //
   // Note that gamma identification code is no longer attributed here: in case a gamma is identified,
   // the IDR->IDcode will be the same general ID code as for all CsI particles,
   // but IDR->IDquality == KVIDGCsI::kICODE10
   //
 
   //perform identification
   KVINDRAIDTelescope::Identify(IDR,x,y);
 
   if(IDR->IDOK)
   {
      if(IDR->Z==0)
      {
         // well-identified particle with Z=0 is a gamma
         IDR->IDquality = KVIDGCsI::kICODE10;
         IDR->Z = 0;
         IDR->A = 0;
         IDR->IDOK = kTRUE;
         IDR->SetComment("gamma");
      }
   }
   else if (IDR->Rejecting_Cut.Contains("gamma")) {
         // gamma identification
         IDR->IDquality = KVIDGCsI::kICODE10;
         IDR->Z = 0;
         IDR->A = 0;
         IDR->IDOK = kTRUE;
         IDR->SetComment("gamma");
      }
      else {
         // rejected by some arbitrary cut line/contour
         IDR->IDquality = KVIDGCsI::kICODE8;
         IDR->Z = 0;
         IDR->A = 0;
         IDR->IDOK = kFALSE;
         IDR->SetComment(Form("rejected by cut:%s", IDR->Rejecting_Cut.Data()));
      }
 
   // set general ID code
   IDR->IDcode = GetIDCode();
 
   return kTRUE;
}
 
 
 
 
void KVIDINDRACsI::SetIdentificationStatus(KVIdentificationResult* IDR, const KVNucleus* nuc)
{
   // For filtering simulations
   //
   // If energy loss in CsI is above threshold for mass identification, we set
   // IDR->Aident=true (and IDR->Zident=true).
   // Otherwise, we just set IDR->Zident=true and A will be given by mass formula for the particle
   //
   // individual thresholds defined for 1H, 2H, 3H, 3He, 4He
   // for A>5 identification if CsI energy > 40 MeV
 
   KVINDRAIDTelescope::SetIdentificationStatus(IDR, nuc);
   if (!IDR->IDOK) return; // check for weird identifications of non-existent nuclei
 
   IDR->Zident = true;
 
   if (IDR->A > 5) {
      if (GetDetector(1)->GetEnergy() > 40)
         IDR->Aident = true;
      return;
   }
   if (fThresMin[IDR->Z - 1][IDR->A - 1] > 0) {
      Bool_t okmass = gRandom->Uniform() < smootherstep(fThresMin[IDR->Z - 1][IDR->A - 1],
                      fThresMax[IDR->Z - 1][IDR->A - 1], GetDetector(1)->GetEnergy());
      if (okmass) {
         IDR->Aident = true;
      }
   }
}
 
 
 
 
 
float KVIDINDRACsI::clamp(float x, float lowerlimit, float upperlimit)
{
   if (x < lowerlimit)
      x = lowerlimit;
   if (x > upperlimit)
      x = upperlimit;
   return x;
}
 
 
 
 
float KVIDINDRACsI::smootherstep(float edge0, float edge1, float x)
{
   // Scale, and clamp x to 0..1 range
   x = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
   // Evaluate polynomial
   return x * x * x * (x * (x * 6 - 15) + 10);
}