00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00033
00034
00044
00045
00046
00047 #include <iostream>
00048 #include <strstream>
00049
00050 #include "TH1D.h"
00051 #include "TTree.h"
00052
00053 #include "globals.hh"
00054 #include "G4DynamicParticle.hh"
00055 #include "G4Event.hh"
00056 #include "G4PrimaryParticle.hh"
00057 #include "G4PrimaryVertex.hh"
00058 #include "G4Step.hh"
00059 #include "G4ThreeVector.hh"
00060 #include "G4Track.hh"
00061 #include "G4TrackStatus.hh"
00062
00063 #include "io/MJLogger.hh"
00064
00065
00066
00067 #include "io/MJOutputSolidBlock.hh"
00068
00069
00070
00071 MJOutputSolidBlock::MJOutputSolidBlock(G4bool isMother):
00072 MJOutputRoot(isMother), fPrimarye(0),
00073 feTrackDistanceHisto(0), feEdepHisto(0),
00074 feVertexDistance(0), feNStepsHisto(0),
00075 feEdepStep(0), feEAtEnd(0), feKEnergyHisto(0), Edep(0)
00076 {
00077 SetSchemaDefined(false);
00078 }
00079
00080
00081
00082 MJOutputSolidBlock::~MJOutputSolidBlock()
00083 {
00084 delete [] Edep;
00085 }
00086
00087
00088
00089 void MJOutputSolidBlock::BeginOfEventAction(const G4Event *event)
00090 {
00091
00092 G4PrimaryParticle *primaryParticle;
00093 G4PrimaryVertex *primaryVertex;
00094 NumberOfParticles = event->GetNumberOfPrimaryVertex();
00095
00096 if(NumberOfParticles > fMaxNumParticles) {
00097 MJLog(warning) << NumberOfParticles << "particles generated. Currently"
00098 << " only " << fMaxNumParticles << " are supported. "
00099 << " EventID : " << event->GetEventID() << endlog;
00100 NumberOfParticles = fMaxNumParticles;
00101 }
00102 MJLog(debugging)<< "Saving " << NumberOfParticles<<" particle(s)." << endlog;
00103
00104
00105 G4int i;
00106 for(i = 0; i < NumberOfParticles; i++) {
00107 primaryVertex = event->GetPrimaryVertex(i);
00108 if(!primaryVertex) {
00109 MJLog(error) << "Bad primary vertex pointer." << endlog;
00110 MJLog(fatal) << endlog;
00111 }
00112 primaryParticle = primaryVertex->GetPrimary();
00113 if(!primaryParticle) {
00114 MJLog(error) << "Bad primary particle pointer." << endlog;
00115 MJLog(fatal) << endlog;
00116 }
00117 Particle[i] = G4ToRoot(primaryParticle->GetPDGcode());
00118 MomentumX[i] =
00119 G4ToRoot(primaryParticle->GetPx() / keV);
00120 MomentumY[i] =
00121 G4ToRoot(primaryParticle->GetPy() / keV);
00122 MomentumZ[i] =
00123 G4ToRoot(primaryParticle->GetPz() / keV);
00124 PositionX[i] = G4ToRoot(primaryVertex->GetX0()/cm);
00125 PositionY[i] = G4ToRoot(primaryVertex->GetY0()/cm);
00126 PositionZ[i] = G4ToRoot(primaryVertex->GetZ0()/cm);
00127 }
00128
00129 primaryVertex = event->GetPrimaryVertex(0);
00130 if(!primaryVertex) {
00131 MJLog(error) << "Bad primary vertex pointer." << endlog;
00132 MJLog(fatal) << endlog;
00133 }
00134 primaryParticle = primaryVertex->GetPrimary();
00135 if(!primaryParticle) {
00136 MJLog(error) << "Bad primary particle pointer." << endlog;
00137 MJLog(fatal) << endlog;
00138 }
00139
00140 if(Edep)
00141 for(i = 0; i < fEdepArraySize; i++)
00142 Edep[i] = 0.0;
00143 else {
00144 MJLog(error) << "Output Array not allocated!" << endlog;
00145 MJLog(fatal) << endlog;
00146 }
00147
00148 EdepOutside = 0.0;
00149 EventNumber = G4ToRoot(event->GetEventID());
00150 fTrackLength = 0.0;
00151 feNSteps = 0;
00152 fPrimarye = 0;
00153 fPrimaryDemiseRecorded = false;
00154 MJLog(debugging) << "BeginOfEventAction finished." << endlog;
00155 }
00156
00157
00158
00159 void MJOutputSolidBlock::BeginOfRunAction()
00160 {
00161 if(IsMother())
00162 OpenRootFile(GetFileName());
00163
00164 feTrackDistanceHisto = new TH1D("h1","Primary e- Track Distance",200,0.,.25);
00165 feEdepHisto = new TH1D("h2","Energy Deposit from Vertex", 200, 0., 0.5);
00166 feVertexDistance = new TH1D("h3","Distance from Vertex", 200, 0., 2.0);
00167 feNStepsHisto = new TH1D("h4","Number of Steps", 400, -0.5, 399.5);
00168 feEdepStep = new TH1D("h5","Edep/step", 200, 0.0, 100.0);
00169 feEAtEnd = new TH1D("h6","E at End", 100, 0.0, 10.0);
00170 feG4Distance = new TH1D("h7","G4 len", 500, 0., 5.0);
00171 feAccTrackLen = new TH1D("h8","Acc Len", 500, 0., 5.0);
00172 feKEnergyHisto = new TH1D("h9","Kinetic Energy", 1000, 0., 1.0);
00173 fEDistance = new TH1D("h10","Distance vs energy loss", 1000, 0., 2.0);
00174
00175 if(!(feTrackDistanceHisto && feEdepHisto && feVertexDistance &&
00176 feNStepsHisto && feEdepStep && feEAtEnd && feG4Distance &&
00177 feAccTrackLen)) {
00178 MJLog(error) << "Could not allocate Histograms! Out of Memory?" << endlog;
00179 MJLog(fatal) << endlog;
00180 }
00181
00182 dX = dY = dZ = 0.001 * cm;
00183 nX = nY = 100;
00184 nZ = 100;
00185 OriginX = OriginY = -0.05 * cm;
00186 OriginZ = -0.002 * cm;
00187 fOrigin = G4ThreeVector(0,0,0);
00188 fMaxNumParticles = 2;
00189 fEdepArraySize = nX * nZ;
00190
00191 Particle = new Int_t[fMaxNumParticles];
00192 MomentumX = new Float_t[fMaxNumParticles];
00193 MomentumY = new Float_t[fMaxNumParticles];
00194 MomentumZ = new Float_t[fMaxNumParticles];
00195 PositionX = new Float_t[fMaxNumParticles];
00196 PositionY = new Float_t[fMaxNumParticles];
00197 PositionZ = new Float_t[fMaxNumParticles];
00198 if(!(Particle && MomentumX && MomentumY && MomentumZ &&
00199 PositionX && PositionY && PositionZ)) {
00200 MJLog(error) << "Could not allocate arrays. Out of Memory?" << endlog;
00201 MJLog(fatal) << endlog;
00202 }
00203
00204 DefineSchema();
00205 MJLog(trace) << "BeginOfRunAction() Finished." << endlog;
00206 }
00207
00208
00209
00210 G4bool MJOutputSolidBlock::ConvertCoordinateToGridPoint(
00211 G4ThreeVector pos, G4int *gridPoint)
00212 {
00213
00214 gridPoint[0] = (G4int)(((Float_t)pos.x() - OriginX) / dX);
00215 gridPoint[1] = (G4int)(((Float_t)pos.y() - OriginY) / dY);
00216 gridPoint[2] = (G4int)(((Float_t)pos.z() - OriginZ) / dZ);
00217
00218
00219 return gridPoint[0] >= 0 && gridPoint[1] >= 0 && gridPoint[2] >= 0 &&
00220 gridPoint[0] < nX && gridPoint[1] < nY && gridPoint[2] < nZ;
00221 }
00222
00223
00224 void MJOutputSolidBlock::DefineSchema()
00225 {
00226 if(!SchemaDefined()){
00227 TTree *nT;
00228 if(fTree)
00229 nT = fTree;
00230 else {
00231 if(!IsMother())
00232 MJLog(warning) << "No tree assigned to child output class." << endlog;
00233 nT = fTree = new TTree("fTree", "LANL Clover Tree");
00234 }
00235
00236
00237 nT->Branch("EventNumber", &EventNumber, "EventNumber/i");
00238 nT->Branch("NumberOfParticles", &NumberOfParticles,"NumberOfParticles/I");
00239
00240 nT->Branch("Particle", Particle, "Particle[NumberOfParticles]/I");
00241
00242 nT->Branch("MomentumX", MomentumX, "MomentumX[NumberOfParticles]/F");
00243 nT->Branch("MomentumY", MomentumY, "MomentumY[NumberOfParticles]/F");
00244 nT->Branch("MomentumZ", MomentumZ, "MomentumZ[NumberOfParticles]/F");
00245
00246 nT->Branch("PositionX", PositionX, "PositionX[NumberOfParticles]/F");
00247 nT->Branch("PositionY", PositionY, "PositionY[NumberOfParticles]/F");
00248 nT->Branch("PositionZ", PositionZ, "PositionZ[NumberOfParticles]/F");
00249
00250
00251
00252 nT->Branch("nX", &nX, "nX/I");
00253 nT->Branch("nY", &nY, "nY/I");
00254 nT->Branch("nZ", &nZ, "nZ/I");
00255
00256
00257 nT->Branch("dX", &dX, "dX/F");
00258 nT->Branch("dY", &dY, "dY/F");
00259 nT->Branch("dZ", &dZ, "dZ/F");
00260
00261
00262 nT->Branch("OriginX", &OriginX, "OriginX/F");
00263 nT->Branch("OriginY", &OriginY, "OriginY/F");
00264 nT->Branch("OriginZ", &OriginZ, "OriginZ/F");
00265
00266
00267 Edep = new Float_t[nX*nZ];
00268 if(!Edep) {
00269 MJLog(error) << "Could not allocate grid array!" << endlog;
00270 MJLog(fatal) << endlog;
00271 }
00272 MJLog(routine) << "Allocated output array of size " << nX*nZ << endlog;
00273
00274
00275
00276
00277 char branchString[20];
00278 strstream branchStream(branchString,sizeof(branchString),
00279 ios::in|ios::out);
00280 branchStream << "Edep[" << nX * nZ << "]/F" << '\0';
00281 nT->Branch("Edep", Edep, branchString);
00282 MJLog(trace) << "Created Branch Edep: " << branchString << endlog;
00283
00284
00285 nT->Branch("EdepOutside", &EdepOutside, "EdepOutside/F");
00286
00287 SetSchemaDefined(true);
00288 }
00289 }
00290
00291
00292
00293
00294 void MJOutputSolidBlock::EndOfEventAction(const G4Event *event)
00295 {
00296 G4int i;
00297 for(i = 0; i < fEdepArraySize; i++)
00298 Edep[i] /= keV;
00299
00300 EdepOutside /= keV;
00301
00302 feNStepsHisto->Fill((Double_t)(G4ToRoot(feNSteps)));
00303 feAccTrackLen->Fill(G4ToRoot(fTrackLength / mm));
00304
00305 if(IsMother())
00306 FillTree();
00307
00308 MJLog(debugging) << "EndOfEventAction Finished." << endlog;
00309 }
00310
00311
00312
00313 void MJOutputSolidBlock::EndOfRunAction()
00314 {
00315 if(IsMother())
00316 CloseRootFile();
00317 MJLog(trace) << "EndOfRunAction finished." << endlog;
00318 }
00319
00320
00321
00322 void MJOutputSolidBlock::RootSteppingAction(const G4Step *step)
00323 {
00324 G4Track *track = step->GetTrack();
00325 G4ThreeVector position = track->GetPosition();
00326 G4int gridPoint[3];
00327 if(ConvertCoordinateToGridPoint(position, gridPoint))
00328 Edep[gridPoint[0] * nZ + gridPoint[2]] += step->GetTotalEnergyDeposit();
00329 else
00330 EdepOutside += step->GetTotalEnergyDeposit();
00331 G4ThreeVector displacement = position - fOrigin;
00332 G4double distancemm = displacement.mag() / mm;
00333 fEDistance->Fill(G4ToRoot(distancemm),
00334 G4ToRoot(step->GetTotalEnergyDeposit() / keV));
00335
00336
00337 if(!fPrimarye)
00338 fPrimarye = track->GetDynamicParticle();
00339 if(track->GetDynamicParticle() == fPrimarye){
00340 if(track->GetTrackStatus() != fAlive && !fPrimaryDemiseRecorded) {
00341 feEAtEnd->Fill(G4ToRoot(fOldKineticEnergy / keV));
00342 feVertexDistance->Fill(G4ToRoot(distancemm));
00343 feG4Distance->Fill(G4ToRoot(track->GetTrackLength() / mm));
00344 fPrimaryDemiseRecorded = true;
00345 } else
00346 fOldKineticEnergy = track->GetKineticEnergy();
00347 feNSteps++;
00348 fTrackLength += step->GetStepLength();
00349 feTrackDistanceHisto->Fill(G4ToRoot(step->GetStepLength() / mm));
00350 feEdepStep->Fill(G4ToRoot(step->GetTotalEnergyDeposit() / keV));
00351 feKEnergyHisto->Fill(G4ToRoot(track->GetKineticEnergy() / keV));
00352 } else {
00353 feEdepHisto->Fill(G4ToRoot(step->GetTotalEnergyDeposit() / keV));
00354 }
00355 }
00356
00357
00358
