HGCal_testbeam/src/PrimaryGeneratorAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 #include "PrimaryGeneratorAction.hh"
0027 #include "PrimaryGeneratorMessenger.hh"
0028
0029 #include "G4Box.hh"
0030 #include "G4LogicalVolume.hh"
0031 #include "G4LogicalVolumeStore.hh"
0032 #include "G4ParticleGun.hh"
0033 #include "G4ParticleTable.hh"
0034 #include "G4RunManager.hh"
0035 #include "Randomize.hh"
0036
0037 #ifdef WITHROOT
0038 #include "G4MTRunManager.hh"
0039 #include "G4Run.hh"
0040 #include "TFile.h"
0041 #include "TTreeReader.h"
0042 #include "TTreeReaderValue.h"
0043 #endif
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 PrimaryGeneratorAction::PrimaryGeneratorAction()
0048     : G4VUserPrimaryGeneratorAction(), fParticleGun(nullptr),
0049       fEnvelopeBox(nullptr) {
0050   G4int n_particle = 1;
0051   fParticleGun = new G4ParticleGun(n_particle);
0052
0053   // default particle kinematic
0054   G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
0055   G4String particleName;
0056   G4ParticleDefinition *particle =
0057       particleTable->FindParticle(particleName = "e+");
0058   fParticleGun->SetParticleDefinition(particle);
0059   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
0060   fParticleGun->SetParticleEnergy(30. * CLHEP::GeV);
0061   fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
0062
0063   fMessenger = new PrimaryGeneratorMessenger(this);
0064 #ifdef WITHROOT
0065   if (fReadInputFile)
0066     OpenInput();
0067 #endif
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 #ifdef WITHROOT
0073 void PrimaryGeneratorAction::OpenInput() {
0074   if (fInputFile != nullptr) {
0075     delete fInputFile;
0076     fInputFile = nullptr;
0077   }
0078   fInputFile = new TFile(fPathInputFile, "READ");
0079   if (fInputFile == nullptr || fInputFile->IsZombie()) {
0080     G4ExceptionDescription msg;
0081     msg << "Input file '" << fPathInputFile << "' cannot be opened.\n";
0082     G4Exception("PrimaryGeneratorAction::GeneratePrimaries()", "WrongInputName",
0083                 FatalErrorInArgument, msg);
0084   }
0085   auto fileDirectory =
0086       dynamic_cast<TDirectory *>(fInputFile->Get("VirtualDetector"));
0087   if (fileDirectory == nullptr) {
0088     G4ExceptionDescription msg;
0089     msg << "Input file '" << fPathInputFile
0090         << "' does not contain 'VirtualDetector' directory.\n";
0091     G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
0092                 "WrongInputDirectory", FatalErrorInArgument, msg);
0093   }
0094   fHgcalReader = new TTreeReader("HGCAL", fileDirectory);
0095   fHgcalEventId = new TTreeReaderValue<Float_t>(*fHgcalReader, "EventID");
0096   fHgcalPdgId = new TTreeReaderValue<Float_t>(*fHgcalReader, "PDGid");
0097   fHgcalPosX = new TTreeReaderValue<Float_t>(*fHgcalReader, "x");
0098   fHgcalPosY = new TTreeReaderValue<Float_t>(*fHgcalReader, "y");
0099   fHgcalMomX = new TTreeReaderValue<Float_t>(*fHgcalReader, "Px");
0100   fHgcalMomY = new TTreeReaderValue<Float_t>(*fHgcalReader, "Py");
0101   fHgcalMomZ = new TTreeReaderValue<Float_t>(*fHgcalReader, "Pz");
0102
0103   return;
0104 }
0105 #endif
0106
0107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0108
0109 PrimaryGeneratorAction::~PrimaryGeneratorAction() {
0110   delete fParticleGun;
0111 #ifdef WITHROOT
0112   delete fInputFile;
0113 #endif
0114 }
0115
0116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0117
0118 void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent) {
0119 #ifdef WITHROOT
0120   if (fReadInputFile) {
0121     // get primary particles from input file
0122     if (fInputFile == nullptr) {
0123       OpenInput();
0124     }
0125     Float_t currentEventId = -1;
0126     while (fHgcalReader->Next()) {
0127       // First particle - set event ID
0128       if (currentEventId == -1) {
0129         currentEventId = **fHgcalEventId;
0130         fEventCounter++;
0131         // check if event should be ignored (is before the start position)
0132         if (fEventCounter < fStartFromEvent)
0133           continue;
0134       } else {
0135         // check if event is to be skipped
0136         if (fEventCounter < fStartFromEvent) {
0137           if (currentEventId != **fHgcalEventId)
0138             currentEventId = -1;
0139           continue;
0140         } else {
0141           // check if current particle belongs to the next event
0142           if (currentEventId != **fHgcalEventId)
0143             break;
0144         }
0145       }
0146       auto vertex = new G4PrimaryVertex(
0147         G4ThreeVector(**fHgcalPosX, **fHgcalPosY, 32.5 * CLHEP::mm), 0 * CLHEP::s);
0148       auto particleDefinition =
0149           G4ParticleTable::GetParticleTable()->FindParticle(**fHgcalPdgId);
0150       auto particle = new G4PrimaryParticle(particleDefinition);
0151       G4ThreeVector momentum(**fHgcalMomX, **fHgcalMomY, **fHgcalMomZ);
0152       particle->SetMomentumDirection(momentum.unit());
0153       particle->SetKineticEnergy(momentum.mag());
0154       vertex->SetPrimary(particle);
0155       anEvent->AddPrimaryVertex(vertex);
0156     }
0157     // if no particles are left - terminate run
0158     if (fEventCounter >= fStartFromEvent &&
0159         anEvent->GetNumberOfPrimaryVertex() == 0) {
0160       G4ExceptionDescription msg;
0161       msg << "Input file does not contain any more events (current event ID = "
0162           << anEvent->GetEventID() << ").\n";
0163       msg << "Run will be aborted with " << anEvent->GetEventID()
0164           << " events processed.\n";
0165       G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
0166                   "EndOfInputFile", JustWarning, msg);
0167       G4RunManager::GetRunManager()->AbortRun();
0168     }
0169     // warn if this is the last event and some events are left in the file
0170     // unprocessed
0171     if (G4RunManager::GetRunManager()
0172             ->GetCurrentRun()
0173             ->GetNumberOfEventToBeProcessed() == anEvent->GetEventID() + 1) {
0174       G4ExceptionDescription msg;
0175       msg << "Input file contains more events than requested in the run.\n";
0176       G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
0177                   "UnusedEventsInInputFile", JustWarning, msg);
0178     }
0179   } else
0180 #endif
0181   {
0182     // Beam position: z
0183     // Use value set by UI command, or set it to edge of the world volume
0184     if (fBeamZ0 == -999 * CLHEP::m) {
0185       G4double worldDZ = 0;
0186       if (!fEnvelopeBox) {
0187         G4LogicalVolume *envLV =
0188             G4LogicalVolumeStore::GetInstance()->GetVolume("World");
0189         if (envLV)
0190           fEnvelopeBox = dynamic_cast<G4Box *>(envLV->GetSolid());
0191       }
0192       if (fEnvelopeBox) {
0193         worldDZ = fEnvelopeBox->GetZHalfLength();
0194         fBeamZ0 = -worldDZ;
0195       } else {
0196         G4ExceptionDescription msg;
0197         msg << "World volume of box shape not found.\n";
0198         msg << "Perhaps you have changed geometry.\n";
0199         msg << "The gun will be place at the center.";
0200         G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
0201                     "WorldVolume", JustWarning, msg);
0202       }
0203     }
0204     // Beam position: x,y
0205     switch (fBeamType) {
0206     case eNone:
0207     default:
0208       fParticleGun->SetParticlePosition(G4ThreeVector(0, 0, fBeamZ0));
0209       break;
0210     case eGaussian:
0211       fParticleGun->SetParticlePosition(
0212           G4ThreeVector(G4RandGauss::shoot(0., fSigmaBeamX),
0213                         G4RandGauss::shoot(0., fSigmaBeamY), fBeamZ0));
0214       break;
0215     case eFlat:
0216       fParticleGun->SetParticlePosition(
0217           G4ThreeVector(G4RandFlat::shoot(-fSigmaBeamX, fSigmaBeamX),
0218                         G4RandFlat::shoot(-fSigmaBeamY, fSigmaBeamY), fBeamZ0));
0219       break;
0220     }
0221     // Particle momentum
0222     if (fMomentumGaussianSpread > 0) {
0223       double energy = fParticleGun->GetParticleEnergy();
0224       energy += G4RandGauss::shoot(0., fMomentumGaussianSpread) * energy;
0225       fParticleGun->SetParticleEnergy(energy);
0226     }
0227     fParticleGun->GeneratePrimaryVertex(anEvent);
0228   }
0229 }
0230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......