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 "Rtypes.h"
00048 #include "TROOT.h"
00049 #include "TFile.h"
00050 #include "TH1.h"
00051 #include "TTree.h"
00052 #include "TObject.h"
00053 #include "TNtuple.h"
00054
00055
00056 #include <strstream>
00057 #include <string>
00058
00059
00060 #include "globals.hh"
00061 #include "G4Run.hh"
00062 #include "G4Event.hh"
00063 #include "G4DynamicParticle.hh"
00064 #include "G4PrimaryParticle.hh"
00065 #include "G4PrimaryVertex.hh"
00066 #include "G4Step.hh"
00067 #include "G4ThreeVector.hh"
00068 #include "G4Track.hh"
00069 #include "G4TrackStatus.hh"
00070
00071
00072 #include "G4SDManager.hh"
00073 #include "G4HCofThisEvent.hh"
00074
00075
00076
00077 #include "gerdageometry/MaGeGeometrySDHit.hh"
00078 #include "gerdaio/MaGeOutputGermaniumArray.hh"
00079 #include "gerdaio/MaGeTrajectory.hh"
00080 #include "io/MJLogger.hh"
00081
00082
00083 MaGeOutputGermaniumArray::MaGeOutputGermaniumArray(G4bool isMother, G4bool needt):
00084 MJOutputRoot(isMother),
00085 NeedTrajectoryInformation(needt)
00086 {
00087 SetSchemaDefined(false);
00088 }
00089
00090
00091 MaGeOutputGermaniumArray::~MaGeOutputGermaniumArray()
00092 {;}
00093
00094 void MaGeOutputGermaniumArray::DefineSchema()
00095 {
00096 if(!SchemaDefined()){
00097 TTree *nT;
00098 if(fTree)
00099 nT = fTree;
00100 else {
00101 if(!IsMother())
00102 MJLog(warning) << "No tree assigned to child output class." << endlog;
00103 nT = fTree = new TTree("fTree", "Germanium Array Tree");
00104 }
00105
00106
00107 nT->Branch("eventnumber",&eventnumber,"eventnumber/I");
00108 nT->Branch("numvertex",&numvertex,"numvertex/I");
00109 nT->Branch("xpos_vertex",xpos_vertex,"xpos_vertex[numvertex]/F");
00110 nT->Branch("ypos_vertex",ypos_vertex,"ypos_vertex[numvertex]/F");
00111 nT->Branch("zpos_vertex",zpos_vertex,"zpos_vertex[numvertex]/F");
00112 nT->Branch("time_vertex",time_vertex,"time_vertex[numvertex]/F");
00113 nT->Branch("numparticle_vertex",numparticle_vertex,"numparticle_vertex[numvertex]/I");
00114 nT->Branch("numparticles", &numparticles, "numparticles/I");
00115 nT->Branch("particle",particle,"particle[numparticles]/I");
00116 nT->Branch("ivertex_particle",ivertex_particle,"ivertex_particle[numparticles]/F");
00117 nT->Branch("px_particle",px_particle,"px_particle[numparticles]/F");
00118 nT->Branch("py_particle",py_particle,"py_particle[numparticles]/F");
00119 nT->Branch("pz_particle",pz_particle,"pz_particle[numparticles]/F");
00120 nT->Branch("pe_particle",pe_particle,"pe_particle[numparticles]/F");
00121 nT->Branch("id_particle",id_particle,"id_particle[numparticles]/I");
00122
00123
00124
00125
00126 nT->Branch("numhits",&numhits,"numhits/I");
00127 nT->Branch("tote_hits",&tote_hits,"tote_hits/F");
00128 nT->Branch("edep_hits",edep_hits,"edep_hits[numhits]/F");
00129 nT->Branch("xpos_hits",xpos_hits,"xpos_hits[numhits]/F");
00130 nT->Branch("ypos_hits",ypos_hits,"ypos_hits[numhits]/F");
00131 nT->Branch("zpos_hits",zpos_hits,"zpos_hits[numhits]/F");
00132 nT->Branch("idge_hits",idge_hits,"idge_hits[numhits]/I");
00133
00134 nT->Branch("numgedet",&numgedet,"numgedet/I");
00135 nT->Branch("idge_gedet",idge_gedet,"idge_gedet[numgedet]/I");
00136 nT->Branch("numhits_gedet",numhits_gedet,"numhits_gedet[numgedet]/I");
00137 nT->Branch("edep_gedet",edep_gedet,"edep_gedet[numgedet]/F");
00138 if (NeedTrajectoryInformation){
00139
00140 nT->Branch("numtrj",&numtrj,"numtrj/I");
00141 nT->Branch("id_trj",id_trj,"id_trj[numtrj]/I");
00142 nT->Branch("pid_trj",pid_trj,"pid_trj[numtrj]/I");
00143 nT->Branch("pdgencode_trj",pdgencode_trj,"pdgencode_trj[numtrj]/I");
00144 nT->Branch("charge_trj",charge_trj,"charge_trj[numtrj]/F");
00145 nT->Branch("px_trj",px_trj,"px_trj[numtrj]/F");
00146 nT->Branch("py_trj",py_trj,"py_trj[numtrj]/F");
00147 nT->Branch("pz_trj",pz_trj,"pz_trj[numtrj]/F");
00148 nT->Branch("npoints_trj",npoints_trj,"npoints_trj[numtrj]/I");
00149 nT->Branch("istart_trj",istart_trj,"istart_trj[numtrj]/I");
00150 nT->Branch("iend_trj",iend_trj,"iend_trj[numtrj]/I");
00151
00152 nT->Branch("numtrjp",&numtrjp,"numtrjp/I");
00153 nT->Branch("x_trjp",x_trjp,"x_trjp[numtrjp]/F");
00154 nT->Branch("y_trjp",y_trjp,"y_trjp[numtrjp]/F");
00155 nT->Branch("z_trjp",z_trjp,"z_trjp[numtrjp]/F");
00156 nT->Branch("de_trjp",de_trjp,"de_trjp[numtrjp]/F");
00157 nT->Branch("steplength_trjp",steplength_trjp,"steplength_trjp[numtrjp]/F");
00158 nT->Branch("insidege_trjp",insidege_trjp,"inside_trjp[numtrjp]/I");
00159 }
00160 SetSchemaDefined(true);
00161 }
00162 }
00163
00164 void MaGeOutputGermaniumArray::BeginOfRunAction()
00165 {
00166 if(IsMother())
00167 OpenRootFile(GetFileName());
00168 DefineSchema();
00169 }
00170
00171
00172 void MaGeOutputGermaniumArray::BeginOfEventAction(const G4Event *event)
00173 {
00174
00175 G4PrimaryParticle *primaryParticle;
00176 G4PrimaryVertex *primaryVertex;
00177 numvertex = event->GetNumberOfPrimaryVertex();
00178
00179 if(numvertex > MAX_NVTX) {
00180 MJLog(warning) << numvertex << "vertex generated. Currently"
00181 << " only " << MAX_NVTX << " are supported. "
00182 << " EventID : " << event->GetEventID() << endlog;
00183 numvertex = MAX_NVTX;
00184 }
00185 MJLog(debugging)<< "Saving " << numvertex<<" vertex(es)." << endlog;
00186
00187 numparticles=0;
00188
00189 for(G4int i = 0; i < numvertex; i++) {
00190 primaryVertex = event->GetPrimaryVertex(i);
00191 if(!primaryVertex) {
00192 MJLog(error) << "Bad primary vertex pointer." << endlog;
00193 MJLog(fatal) << endlog;
00194 }
00195 xpos_vertex[i] = G4ToRoot(primaryVertex->GetX0()/cm);
00196 ypos_vertex[i] = G4ToRoot(primaryVertex->GetY0()/cm);
00197 zpos_vertex[i] = G4ToRoot(primaryVertex->GetZ0()/cm);
00198 time_vertex[i] = G4ToRoot(primaryVertex->GetT0()/ns);
00199 numparticle_vertex[i] = primaryVertex->GetNumberOfParticle();
00200
00201 for (G4int j=0; j<numparticle_vertex[i]; j++) {
00202 primaryParticle = primaryVertex->GetPrimary(j);
00203 if(!primaryParticle) {
00204 MJLog(error) << "Bad primary particle pointer." << endlog;
00205 MJLog(fatal) << endlog;
00206 }
00207 particle[numparticles] = G4ToRoot(primaryParticle->GetPDGcode());
00208 px_particle[numparticles] =
00209 G4ToRoot(primaryParticle->GetPx() / keV);
00210 py_particle[numparticles] =
00211 G4ToRoot(primaryParticle->GetPy() / keV);
00212 pz_particle[numparticles] =
00213 G4ToRoot(primaryParticle->GetPz() / keV);
00214 ivertex_particle[numparticles]=i;
00215 numparticles++;
00216 }
00217 }
00218 eventnumber = G4ToRoot(event->GetEventID());
00219
00220 }
00221
00222 void MaGeOutputGermaniumArray::EndOfEventAction(const G4Event *evt)
00223 {
00224
00225
00226
00227
00228
00229 G4SDManager * SDman = G4SDManager::GetSDMpointer();
00230 G4int geCollID = SDman->GetCollectionID("HitsCollection");
00231
00232
00233 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00234 MaGeGeometrySDHitsCollection* GCHC = 0;
00235 if (HCE) GCHC = (MaGeGeometrySDHitsCollection*)(HCE->GetHC(geCollID));
00236
00237 G4ThreeVector tempvector;
00238
00239 if (GCHC)
00240 {
00241 tote_hits = 0.0;
00242 for (G4int i=0; i<MAX_NGEDET; i++) {
00243 edep_allgedet[i]=0.0;
00244 numhits_allgedet[i]=0;
00245 }
00246 numhits = GCHC->entries();
00247 for (G4int i=0; i<numhits; i++) {
00248 tote_hits += (*GCHC)[i]->GetEdep();
00249 edep_hits[i] = (*GCHC)[i]->GetEdep();
00250 tempvector = (*GCHC)[i]->GetPos();
00251 xpos_hits[i] = tempvector.getX();
00252 ypos_hits[i] = tempvector.getY();
00253 zpos_hits[i] = tempvector.getZ();
00254 idge_hits[i] = (*GCHC)[i]->GetCopynumber();
00255 edep_allgedet[idge_hits[i]] += edep_hits[i];
00256 numhits_allgedet[idge_hits[i]]++;
00257
00258
00259
00260 }
00261 }
00262
00263 numgedet=0;
00264 for (G4int i=0; i<MAX_NGEDET; i++) {
00265 if (edep_allgedet[i]>0.0) {
00266 idge_gedet[numgedet]=i;
00267 edep_gedet[numgedet]=edep_allgedet[i];
00268 numhits_gedet[numgedet]=numhits_allgedet[i];
00269 numgedet++;
00270 }
00271 }
00272
00273
00274
00275
00276 if (NeedTrajectoryInformation){
00277 G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer();
00278 numtrj = 0;
00279 numtrjp = 0;
00280 G4String volnametemp;
00281 if (trajectoryContainer) numtrj = trajectoryContainer->entries();
00282 for (G4int i=0; i<numtrj; i++) {
00283 MaGeTrajectory* trj = (MaGeTrajectory*)
00284 ((*(evt->GetTrajectoryContainer()))[i]);
00285 id_trj[i] =trj->GetTrackID();
00286 pid_trj[i] =trj->GetParentID();
00287 pdgencode_trj[i] =trj->GetPDGEncoding();
00288 charge_trj[i] =trj->GetCharge();
00289 px_trj[i] =trj->GetInitialMomentum().x();
00290 py_trj[i] =trj->GetInitialMomentum().y();
00291 pz_trj[i] =trj->GetInitialMomentum().z();
00292 npoints_trj[i] =trj->GetPointEntries();
00293 istart_trj[i] =numtrjp;
00294 numtrjp +=npoints_trj[i];
00295 iend_trj[i] =numtrjp-1;
00296 for (G4int j=0; j<npoints_trj[i]; j++) {
00297 MaGeTrajectoryPoint* trjp = (MaGeTrajectoryPoint*)(trj->GetPoint(j));
00298 G4int jj=j+istart_trj[i];
00299 x_trjp[jj] =trjp->GetPosition().x();
00300 y_trjp[jj] =trjp->GetPosition().y();
00301 z_trjp[jj] =trjp->GetPosition().z();
00302 de_trjp[jj] =trjp->GetEnergyLost();
00303 steplength_trjp[jj]=trjp->GetStepLength();
00304 volnametemp =trjp->GetVolumeName();
00305 insidege_trjp[jj]=0;
00306 if (volnametemp.find("Ge")!=string::npos) {insidege_trjp[jj]=1;}
00307 if (j==0) {insidege_trjp[jj]=1;}
00308
00309 }
00310 }
00311 }
00312
00313 if (IsMother()) FillTree();
00314 MJLog(debugging) << "EndOfEventAction Finished." << endlog;
00315
00316 }
00317
00318 void MaGeOutputGermaniumArray::EndOfRunAction()
00319 {
00320 if(IsMother())
00321 CloseRootFile();
00322 MJLog(trace) << "EndOfRunAction finished." << endlog;
00323 }
00324
00325 void MaGeOutputGermaniumArray::RootSteppingAction(const G4Step *step)
00326 {;}
00327
