~chris-rogers/maus/1312

« back to all changes in this revision

Viewing changes to src/common_cpp/Simulation/MAUSPhysicsList.cc

• browse files at revision 659.1.93
• compare with another revision
• download diff
• download tarball
• view history from revision 659.1.93

Show diffs side-by-side

added

removed

src/common_cpp/Simulation/MAUSPhysicsList.cc

23

23

 

24

24

#include "json/json.h"

25

25

 

26

 

 

27

 

 

28

26

#include "Geant4/G4ProcessManager.hh"

29

27

#include "Geant4/G4ProcessTable.hh"

30

28

#include "Geant4/G4ProcessVector.hh"

42

40

#include "Geant4/globals.hh"

43

41

#include "Geant4/G4StepLimiter.hh"

44

42

#include "Geant4/G4UserSpecialCuts.hh"

 

43

#include "Geant4/G4ProductionCuts.hh"

45

44

#include "Geant4/G4UImanager.hh"

46

45

#include "Geant4/G4PhysListFactory.hh"

47

46

 

 

47

#include "Geant4/G4Region.hh"

 

48

#include "Geant4/G4RegionStore.hh"

 

49

 

48

50

#include "Interface/Squeak.hh"

 

51

#include "Interface/STLUtils.hh"

49

52

 

50

53

#include "src/common_cpp/Utils/Globals.hh"

51

54

#include "src/common_cpp/Utils/JsonWrapper.hh"

52

55

#include "src/common_cpp/Simulation/MAUSPhysicsList.hh"

 

56

#include "src/common_cpp/Simulation/DetectorConstruction.hh"

 

57

#include "src/common_cpp/Simulation/MAUSGeant4Manager.hh"

53

58

 

54

59

namespace MAUS {

55

60

 

155

160

 

156

161

 

157

162

void MAUSPhysicsList::SetEnergyLoss(eloss eLossModel) {

158

 

  double       cutDouble = _productionThreshold;

159

 

  std::stringstream cutStream;

 

163

  double cutDouble = _productionThreshold;

160

164

  std::string  elossActive = "activate";

161

165

  std::string  flucActive  = "true";

162

166

  switch (eLossModel) {

176

180

      cutDouble   = 1e12*mm;

177

181

      break;

178

182

  }

179

 

  cutStream << cutDouble;

180

183

  std::vector<std::string> uiCommand;

181

 

  uiCommand.push_back("/run/setCut "+cutStream.str());

 

184

  if (eLossModel != energyStraggling)

 

185

    uiCommand.push_back("/run/setCut "+STLUtils::ToString(cutDouble));

182

186

  uiCommand.push_back("/process/eLoss/fluct "+flucActive);

183

187

  for (int i = 0; i < _nELossNames; i++)

184

188

    uiCommand.push_back("/process/"+elossActive+" "+_eLossNames[i]);

186

190

  for (size_t i = 0; i < uiCommand.size(); i++) {

187

191

    UIApplyCommand(uiCommand[i]);

188

192

  }

 

193

 

 

194

  // set by volume

 

195

  MAUSGeant4Manager* g4man = MAUSGeant4Manager::GetInstance();

 

196

  std::vector<std::string> regions = g4man->GetGeometry()->GetRegions();

 

197

  for (size_t i = 0; i < regions.size() && eLossModel == energyStraggling; ++i) {

 

198

      std::map<std::string, double> thresholds;

 

199

      if (_fineGrainedProductionThreshold.find(regions[i]) !=

 

200

                                        _fineGrainedProductionThreshold.end()) {

 

201

          thresholds = _fineGrainedProductionThreshold[regions[i]];

 

202

      }

 

203

      SetProductionThresholdByVolume(regions[i], cutDouble, thresholds);

 

204

  }

189

205

}

190

206

 

191

207

 

368

384

                                          JsonWrapper::realValue).asDouble();

369

385

    _productionThreshold = JsonWrapper::GetProperty(dc, "production_threshold",

370

386

                                          JsonWrapper::realValue).asDouble();

 

387

    Json::Value jsonFineGrained = JsonWrapper::GetProperty(dc,

 

388

                                          "fine_grained_production_threshold",

 

389

                                          JsonWrapper::objectValue);

 

390

    // loop over json {region1:{pid1:cut1, pid2:cut2}, region2:{...}, ...}

 

391

    std::vector<std::string> names = jsonFineGrained.getMemberNames();

 

392

    for (size_t i = 0; i < names.size(); ++i) {

 

393

        Json::Value jsonCuts = JsonWrapper::GetProperty(

 

394

                              jsonFineGrained,

 

395

                              names[i],

 

396

                              JsonWrapper::objectValue);

 

397

        std::vector<std::string> pidStrings = jsonCuts.getMemberNames();

 

398

        std::map<std::string, double> cuts;

 

399

        for (size_t j = 0; j < pidStrings.size(); ++j) {

 

400

            cuts[pidStrings[j]] = JsonWrapper::GetProperty(jsonCuts,

 

401

                                          pidStrings[j],

 

402

                                          JsonWrapper::realValue).asDouble();

 

403

        }

 

404

        _fineGrainedProductionThreshold[names[i]] = cuts;

 

405

    }

371

406

}

372

407

 

373

408

void MAUSPhysicsList::UIApplyCommand(std::string command) {

375

410

        << "Apply G4UI command: " << command << std::endl;

376

411

    G4UImanager::GetUIpointer()->ApplyCommand(command);

377

412

}

 

413

 

 

414

void MAUSPhysicsList::SetProductionThresholdByVolume(

 

415

                          std::string volumeName,

 

416

                          double defaultProductionThreshold,

 

417

                          std::map<std::string, double> particleIdToThreshold) {

 

418

    G4Region* region = G4RegionStore::GetInstance()->GetRegion(volumeName);

 

419

    if (region == NULL) {

 

420

        throw MAUS::Exception(Exception::recoverable,

 

421

                              "Failed to find region "+volumeName+" for G4Cuts",

 

422

                              "MAUSPhysicsList::SetProductionThresholdByVolume");

 

423

    }

 

424

    G4ProductionCuts* cuts = new G4ProductionCuts();

 

425

    Squeak::mout(Squeak::debug) << "Cuts for Volume " << volumeName << std::endl;

 

426

    G4ParticleTable* ptable = G4ParticleTable::GetParticleTable();

 

427

    if (particleIdToThreshold.find("default") != particleIdToThreshold.end()) {

 

428

        defaultProductionThreshold = particleIdToThreshold["default"];

 

429

    }

 

430

 

 

431

    for (int i = 0; i < ptable->size(); ++i) {

 

432

        int pid = ptable->GetParticle(i)->GetPDGEncoding();

 

433

        std::string nameString = ptable->GetParticle(i)->GetParticleName();

 

434

        std::string pidString = STLUtils::ToString(pid);

 

435

        double cut = 0.;

 

436

        if (particleIdToThreshold.find(pidString) != particleIdToThreshold.end()) {

 

437

            cut = particleIdToThreshold[pidString];

 

438

        } else if (particleIdToThreshold.find(nameString) != particleIdToThreshold.end()) {

 

439

            cut = particleIdToThreshold[nameString];

 

440

        } else {

 

441

            cut = defaultProductionThreshold;

 

442

        }

 

443

        if (cut > 0.) {

 

444

            Squeak::mout(Squeak::debug) << "    " << pid << " " << cut << std::endl;

 

445

            cuts->SetProductionCut(cut, pid);

 

446

        } else {

 

447

            Squeak::mout(Squeak::debug) << "    " << pid

 

448

                                        << " -ve so not set" << std::endl;

 

449

        }

 

450

    }

 

451

    region->SetProductionCuts(cuts);

 

452

}

378

453

}

379

454

 

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1