lAr_calorimeter/src/FCALPrimaryGeneratorAction.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 //
0027 //
0028
0029 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0031
0032 #include <fstream>
0033 #include <cstdlib>
0034
0035 #include "FCALPrimaryGeneratorAction.hh"
0036
0037 #include "G4SystemOfUnits.hh"
0038 #include "G4Event.hh"
0039 #include "G4ParticleGun.hh"
0040 #include "G4ParticleTable.hh"
0041 #include "G4ParticleDefinition.hh"
0042 #include "Randomize.hh"
0043 #include "G4DataVector.hh"
0044 #include "G4AutoLock.hh"
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0047
0048
0049 // Migration to MT: there is a single input file that is read by all threads.
0050 // The idea is that the events are read by a single thread and processed
0051 // by all threads. Threads ask for the next ID to be processed. When
0052 // events are all processed we start over from the beginning of the file
0053 namespace {
0054     G4bool isFileRead = false;
0055     G4Mutex mFileRead = G4MUTEX_INITIALIZER;
0056     //Primary kinematics
0057     G4DataVector fX;
0058     G4DataVector fY;
0059     G4DataVector fZ;
0060     G4DataVector fCosX;
0061     G4DataVector fCosY;
0062     G4DataVector fCosZ;
0063     size_t nextEventId = 0;
0064     G4Mutex mNextEventId = G4MUTEX_INITIALIZER;
0065
0066     size_t GetNextId() {
0067         G4AutoLock l(&mNextEventId);
0068         if ( nextEventId >= fX.size() ) //file data are over,  restart file
0069             {
0070                 G4Exception("FCALPrimaryGeneratorAction::GeneratePrimaries","lAr002",
0071                             JustWarning,"Data file with kinematics is over, restart it");
0072                 nextEventId=0;
0073             }
0074         return nextEventId++;
0075     }
0076
0077     void ReadKinematicFromFile(G4double energy) {
0078         //Only one thread shoud read input file
0079         G4AutoLock l(&mFileRead);
0080         if ( isFileRead ) return;
0081         // Read Kinematics from file
0082         G4String file_name = "data-tracks/tracks-80GeV.dat";
0083         if (energy < 30*GeV)
0084             file_name = "data-tracks/tracks-20GeV.dat";
0085         else if (energy < 50*GeV)
0086             file_name = "data-tracks/tracks-40GeV.dat";
0087         else if (energy < 70*GeV)
0088             file_name = "data-tracks/tracks-60GeV.dat";
0089         else if (energy < 90*GeV)
0090             file_name = "data-tracks/tracks-80GeV.dat";
0091         else if (energy < 150*GeV)
0092             file_name = "data-tracks/tracks-120GeV.dat";
0093         else
0094             file_name = "data-tracks/tracks-200GeV.dat";
0095         std::ifstream Traks_file(file_name);
0096         if(!Traks_file)
0097         {
0098             G4ExceptionDescription ed;
0099             ed << "Failed to open file " << file_name << G4endl;
0100             G4Exception("FCALPrimaryGeneratorAction::FCALPrimaryGeneratorAction()",
0101                         "lAr001",FatalException,ed);
0102         }
0103         G4double xx=0,yy=0,zz=0,c1=0,c2=0,c3=0;
0104         G4int iev = 0;
0105         while(!(Traks_file.eof())) {
0106             Traks_file >> iev >> xx >> yy >> zz >> c1 >> c2 >> c3;
0107             fX.push_back(xx*cm);
0108             fY.push_back(yy*cm);
0109             fZ.push_back(zz*cm);
0110             fCosX.push_back(c1);
0111             fCosY.push_back(c2);
0112             fCosZ.push_back(c3);
0113         }
0114         G4cout << "Read " << fX.size() << " events from file " << file_name << G4endl;
0115         isFileRead= true;
0116         Traks_file.close();
0117         return;
0118     }
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0122
0123 FCALPrimaryGeneratorAction::FCALPrimaryGeneratorAction() :
0124   fVerbosity(0)
0125 {
0126   particleGun  = new G4ParticleGun();
0127
0128   // default Particle
0129   G4String particleName;
0130   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0131   G4ParticleDefinition* particle = particleTable->FindParticle(particleName="e-");
0132   particleGun->SetParticleDefinition(particle);
0133
0134   // default Energy
0135   particleGun->SetParticleEnergy(20*GeV);
0136 }
0137
0138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0139
0140 FCALPrimaryGeneratorAction::~FCALPrimaryGeneratorAction()
0141 {
0142   delete particleGun;
0143 }
0144
0145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0146
0147 void FCALPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0148 {
0149   //this function is called at the begining of event
0150   ReadKinematicFromFile(particleGun->GetParticleEnergy());
0151
0152     //Get next event to be processed
0153     size_t nEvent = GetNextId();
0154     particleGun->SetParticlePosition(G4ThreeVector(fX[nEvent],fY[nEvent],fZ[nEvent]));
0155   particleGun->SetParticleMomentumDirection(G4ThreeVector(-1.0*fCosX[nEvent],
0156                               fCosY[nEvent],
0157                               -1.0*fCosZ[nEvent]));
0158
0159   particleGun->GeneratePrimaryVertex(anEvent);
0160
0161   if (fVerbosity)
0162     {
0163         G4cout<< " Event  "<<anEvent->GetEventID()<< " Generated Vertex : "
0164             <<anEvent->GetEventID() <<" (x,y,z)=(" << fX[nEvent] << ","
0165             <<fY[nEvent] << "," << fZ[nEvent]<< ") (cosX,cosY,cosZ)=("
0166             << -1.*fCosX[nEvent] << "," << fCosY[nEvent]
0167             <<"," << -1.*fCosZ[nEvent] << ")"<<G4endl;
0168     }
0169
0170 }
0171
0172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0173
0174