TGeant4/TG4G3CutVector.cxx

   1 // $Id$
   2 // Category: global
   3 //
   4 // Author: I. Hrivnacova
   5 //
   6 // Class TG4G3CutVector
   7 // --------------------
   8 // See the class description in the header file.
   9
  10 #include "TG4G3CutVector.h"
  11 #include "TG4G3Defaults.h"
  12
  13 #include <G4Track.hh>
  14 #include <G4ParticleDefinition.hh>
  15 #include <G4VProcess.hh>
  16
  17 #include <g4std/strstream>
  18
  19 const G4double  TG4G3CutVector::fgkDCUTEOff  = 10. * TeV;
  20 const G4double  TG4G3CutVector::fgkDCUTMOff  = 10. * TeV;
  21 const G4double  TG4G3CutVector::fgkTolerance =  1. * keV;
  22 TG4StringVector TG4G3CutVector::fgCutNameVector;
  23
  24 //_____________________________________________________________________________
  25 TG4G3CutVector::TG4G3CutVector()
  26   : fDeltaRaysOn(true)
  27 {
  28   // fill name vector
  29   if (fgCutNameVector.size() == 0) FillCutNameVector();
  30
  31   // initialize fCutVector
  32   for (G4int i=0; i<kTOFMAX; i++) fCutVector.push_back(0.);
  33   fCutVector.push_back(DBL_MAX);
  34 }
  35
  36 //_____________________________________________________________________________
  37 TG4G3CutVector::TG4G3CutVector(const TG4G3CutVector& right)
  38   : fCutVector(right.fCutVector.size())
  39 {
  40   // copy stuff
  41   *this = right;
  42 }
  43
  44 //_____________________________________________________________________________
  45 TG4G3CutVector::~TG4G3CutVector() {
  46 //
  47 }
  48
  49 // operators
  50
  51 //_____________________________________________________________________________
  52 TG4G3CutVector& TG4G3CutVector::operator=(const TG4G3CutVector& right)
  53 {
  54   // check assignement to self
  55   if (this == &right) return *this;
  56
  57   // copy fCutVector
  58   for (G4int i=0; i<kNoG3Cuts; i++) fCutVector[i] = right.fCutVector[i];
  59
  60   fDeltaRaysOn = right.fDeltaRaysOn;
  61
  62   return *this;
  63 }
  64
  65 //_____________________________________________________________________________
  66 G4double TG4G3CutVector::operator[](G4int index) const
  67 {
  68 //
  69   if (index < kNoG3Cuts)
  70     return fCutVector[index];
  71   else {
  72     TG4Globals::Exception(
  73       "TG4G3CutVector::operator[]: index out of the vector scope");
  74     return 0.;
  75   }
  76 }
  77
  78 // private methods
  79
  80 //_____________________________________________________________________________
  81 void TG4G3CutVector::FillCutNameVector()
  82 {
  83 // Defines fCutNameVector.
  84 // ---
  85
  86   fgCutNameVector.push_back("CUTGAM");
  87   fgCutNameVector.push_back("CUTELE");
  88   fgCutNameVector.push_back("CUTNEU");
  89   fgCutNameVector.push_back("CUTHAD");
  90   fgCutNameVector.push_back("CUTMUO");
  91   fgCutNameVector.push_back("BCUTE");
  92   fgCutNameVector.push_back("BCUTM");
  93   fgCutNameVector.push_back("DCUTE");
  94   fgCutNameVector.push_back("DCUTM");
  95   fgCutNameVector.push_back("PPCUTM");
  96   fgCutNameVector.push_back("TOFMAX");
  97   fgCutNameVector.push_back("NONE");
  98 }
  99
 100 // public methods
 101
 102 //_____________________________________________________________________________
 103 TG4G3Cut TG4G3CutVector::GetCut(const G4String& cutName)
 104 {
 105 // Retrieves corresponding TG4G3Cut constant from the cutName.
 106 // ---
 107
 108   if      (cutName == fgCutNameVector[kCUTGAM]) return kCUTGAM;
 109   else if (cutName == fgCutNameVector[kCUTELE]) return kCUTELE;
 110   else if (cutName == fgCutNameVector[kCUTNEU]) return kCUTNEU;
 111   else if (cutName == fgCutNameVector[kCUTHAD]) return kCUTHAD;
 112   else if (cutName == fgCutNameVector[kCUTMUO]) return kCUTMUO;
 113   else if (cutName == fgCutNameVector[kBCUTE])  return kBCUTE;
 114   else if (cutName == fgCutNameVector[kBCUTM])  return kBCUTM;
 115   else if (cutName == fgCutNameVector[kDCUTE])  return kDCUTE;
 116   else if (cutName == fgCutNameVector[kDCUTM])  return kDCUTM;
 117   else if (cutName == fgCutNameVector[kPPCUTM]) return kPPCUTM;
 118   else if (cutName == fgCutNameVector[kTOFMAX]) return kTOFMAX;
 119   else return kNoG3Cuts;
 120 }
 121
 122 //_____________________________________________________________________________
 123 const G4String& TG4G3CutVector::GetCutName(TG4G3Cut cut)
 124 {
 125 // Returns name of a specified cut.
 126 // ---
 127
 128   // fill name vector
 129   if (fgCutNameVector.size() == 0) TG4G3CutVector::FillCutNameVector();
 130
 131   return fgCutNameVector[cut];
 132 }
 133
 134 //_____________________________________________________________________________
 135 void TG4G3CutVector::SetCut(TG4G3Cut cut, G4double cutValue)
 136 {
 137 // Sets the cutValue for the specified cut.
 138 // ---
 139
 140   if (cut>=kNoG3Cuts) {
 141     TG4Globals::Exception(
 142       "TG4G3CutVector::SetG3Cut: Inconsistent cut.");
 143   }
 144
 145   fCutVector[cut] = cutValue;
 146 }
 147
 148 //_____________________________________________________________________________
 149 void TG4G3CutVector::SetG3Defaults()
 150 {
 151 // Sets G3 default values for all cuts.
 152 // ---
 153
 154   for (G4int i=0; i<kNoG3Cuts; i++)
 155     fCutVector[i] = TG4G3Defaults::Instance()->CutValue(i);
 156 }
 157
 158 //_____________________________________________________________________________
 159 G4bool TG4G3CutVector::Update(const TG4G3CutVector& vector)
 160 {
 161 // Update only those values that have not yet been set ( =0.)
 162 // with given vector.
 163 // Returns true if some value was modified.
 164 // ---
 165
 166   G4bool result = false;
 167
 168   for (G4int i=0; i<kNoG3Cuts; i++)
 169     if (fCutVector[i] < fgkTolerance ) {
 170        fCutVector[i] = vector[i];
 171        result = true;
 172     }
 173
 174   return result;
 175 }
 176
 177 //_____________________________________________________________________________
 178 G4String TG4G3CutVector::Format() const
 179 {
 180 // Formats the output into a string.
 181 // ---
 182
 183   strstream tmpStream;
 184   tmpStream << "  G3 cut vector:" << G4endl;
 185   if (fDeltaRaysOn)
 186     tmpStream << "    Delta rays On" << G4endl;
 187   else
 188     tmpStream << "    Delta rays Off" << G4endl;
 189
 190   // energy cuts
 191   for (G4int i=0; i<kTOFMAX; i++) {
 192
 193     tmpStream << "    " << fgCutNameVector[i] << " cut value (MeV): ";
 194
 195     if (i == kDCUTE && !fDeltaRaysOn)
 196       tmpStream << fgkDCUTEOff/MeV << G4endl;
 197     else if (i == kDCUTM && !fDeltaRaysOn)
 198       tmpStream << fgkDCUTMOff/MeV << G4endl;
 199     else
 200       tmpStream << fCutVector[i]/MeV << G4endl;
 201   }
 202
 203   // time cut
 204   tmpStream << "    " << fgCutNameVector[kTOFMAX] << " cut value (s):   "
 205             << fCutVector[kTOFMAX]/s << G4endl;
 206
 207   return tmpStream.str();
 208 }
 209
 210 //_____________________________________________________________________________
 211 void TG4G3CutVector::Print() const
 212 {
 213 // Prints the cuts.
 214 // ---
 215
 216   G4cout << Format();
 217 }
 218
 219 //_____________________________________________________________________________
 220 G4double TG4G3CutVector::GetMinEkineForGamma(const G4Track& track) const
 221 {
 222 // Returns the cut value for gamma.
 223 // (Cut is not applied for "opticalphoton"
 224 //  as it is treated in G4 as a particle different
 225 //  from "gamma" in G4.)
 226 // ---
 227
 228   const G4VProcess* kpCreatorProcess = track.GetCreatorProcess();
 229   G4String processName = "";
 230   if (kpCreatorProcess) processName = kpCreatorProcess->GetProcessName();
 231
 232   if ((processName == "eBrem") || (processName == "IeBrem")) {
 233     return fCutVector[kBCUTE];
 234   }
 235   else if ((processName == "MuBrems") || (processName == "IMuBrems")) {
 236            //(processName == "//hBrems")|| (processName == "//IhBrems")
 237            // hadron Brehmstrahlung is not defined in G4
 238     return fCutVector[kBCUTM];
 239   }
 240   else {
 241     return fCutVector[kCUTGAM];
 242   }
 243 }
 244
 245 //_____________________________________________________________________________
 246 G4double TG4G3CutVector::GetMinEkineForElectron(const G4Track& track) const
 247 {
 248 // Returns the cut value for e-.
 249 // ---
 250
 251   const G4VProcess* kpCreatorProcess = track.GetCreatorProcess();
 252   G4String processName = "";
 253   if (kpCreatorProcess) processName = kpCreatorProcess->GetProcessName();
 254
 255   if ((processName == "eIoni") || (processName == "IeIoni")) {
 256     // delta rays by e-, e+
 257
 258     if (fDeltaRaysOn)
 259       return fCutVector[kDCUTE];
 260     else
 261       return fgkDCUTEOff;
 262   }
 263   else if ((processName == "MuIoni") || (processName == "IMuIoni") ||
 264            (processName == "hIoni")  || (processName == "IhIoni")) {
 265     // delta rays by other particles (mu, hadron)
 266
 267     if (fDeltaRaysOn)
 268       return fCutVector[kDCUTM];
 269     else
 270       return fgkDCUTMOff;
 271   }
 272   else if (processName == "MuPairProd") {
 273     //direct pair production by muons
 274
 275     return fCutVector[kPPCUTM];
 276   }
 277   else {
 278     return fCutVector[kCUTELE];
 279   }
 280 }
 281
 282 //_____________________________________________________________________________
 283 G4double TG4G3CutVector::GetMinEkineForEplus(const G4Track& track) const
 284 {
 285 // Returns the cut value for e+.
 286 // ---
 287
 288   const G4VProcess* kpCreatorProcess = track.GetCreatorProcess();
 289   G4String processName = "";
 290   if (kpCreatorProcess) processName = kpCreatorProcess->GetProcessName();
 291
 292   if (processName == "MuPairProd") {
 293     //direct pair production by muons
 294     return fCutVector[kPPCUTM];
 295   }
 296   else
 297     return fCutVector[kCUTELE];
 298 }
 299
 300 //_____________________________________________________________________________
 301 G4double TG4G3CutVector::GetMinEkineForChargedHadron(const G4Track& track) const
 302 {
 303 // Returns the cut value for charged hadron.
 304 // ---
 305
 306   return fCutVector[kCUTHAD];
 307 }
 308
 309 //_____________________________________________________________________________
 310 G4double TG4G3CutVector::GetMinEkineForNeutralHadron(const G4Track& track) const
 311 {
 312 // Returns the cut value for neutral hadron.
 313 // ---
 314
 315   return fCutVector[kCUTNEU];
 316 }
 317
 318 //_____________________________________________________________________________
 319 G4double TG4G3CutVector::GetMinEkineForMuon(const G4Track& track) const
 320 {
 321 // Returns the cut value for neutral muon.
 322 // ---
 323
 324   return fCutVector[kCUTMUO];
 325 }
 326
 327 //_____________________________________________________________________________
 328 G4double TG4G3CutVector::GetMinEkineForOther(const G4Track& track) const
 329 {
 330 // Returns 0.
 331 // ---
 332
 333   return 0.;
 334 }
 335 //_____________________________________________________________________________
 336 G4bool TG4G3CutVector::IsCut() const
 337 {
 338 // Returns true if any of cuts is set.
 339 // ---
 340
 341   for (G4int i=0; i<kNoG3Cuts; i++)
 342     if (fCutVector[i] > fgkTolerance )  return true;
 343
 344   return false;
 345 }