[286] | 1 | #include "ELYSE/EventAction.hh" |
---|
| 2 | |
---|
| 3 | //Geant4 |
---|
| 4 | #include "G4Event.hh" |
---|
| 5 | #include "G4RunManager.hh" |
---|
| 6 | #include "G4EventManager.hh" |
---|
| 7 | #include "G4UImanager.hh" |
---|
| 8 | #include "G4TrajectoryContainer.hh" |
---|
| 9 | #include "G4VVisManager.hh" |
---|
| 10 | #include "G4ios.hh" |
---|
| 11 | #include "globals.hh" |
---|
| 12 | #include "G4ThreeVector.hh" |
---|
| 13 | #include "G4TransportationManager.hh" |
---|
| 14 | #include "G4Navigator.hh" |
---|
| 15 | #include "G4SDManager.hh" |
---|
| 16 | #include "G4DigiManager.hh" |
---|
| 17 | #include "G4UnitsTable.hh" |
---|
| 18 | #include "G4UIcmdWith3VectorAndUnit.hh" |
---|
| 19 | |
---|
| 20 | //std |
---|
| 21 | #include <set> |
---|
| 22 | #include <iomanip> |
---|
| 23 | #include <string> |
---|
| 24 | #include <vector> |
---|
| 25 | #include <stdlib.h> //abs |
---|
| 26 | |
---|
| 27 | // AIDA : |
---|
| 28 | #include <AIDA/ITree.h> |
---|
| 29 | #include <AIDA/IManagedObject.h> |
---|
| 30 | #include <AIDA/ITuple.h> |
---|
| 31 | |
---|
| 32 | //ELYSE |
---|
| 33 | #include "ELYSE/Cast.hh" |
---|
| 34 | #include "ELYSE/Analysis.hh" |
---|
| 35 | #include "ELYSE/Trajectory.hh" |
---|
| 36 | #include "ELYSE/RunAction.hh" |
---|
| 37 | #include "ELYSE/PrimaryGeneratorAction.hh" |
---|
[294] | 38 | //#include "ELYSE/TrappVolHit.hh" |
---|
[286] | 39 | #include "ELYSE/DetectorConstruction.hh" |
---|
| 40 | |
---|
| 41 | |
---|
| 42 | ELYSE::EventAction::EventAction(ELYSE::Analysis& aAnalysis, |
---|
| 43 | ELYSE::RunAction* myRun, |
---|
| 44 | ELYSE::DetectorConstruction* myDetector, |
---|
| 45 | ELYSE::PrimaryGeneratorAction* myGenerator) |
---|
| 46 | :fAnalysis(aAnalysis), |
---|
| 47 | runAction(myRun), |
---|
| 48 | generatorAction(myGenerator), |
---|
| 49 | detectorConstructor(myDetector), |
---|
| 50 | eventTuple(0) { |
---|
| 51 | |
---|
| 52 | //Get User Histo pointers |
---|
| 53 | AIDA::ITree* usrTree = aAnalysis.tree(); //get the tree |
---|
| 54 | if (!usrTree) { |
---|
| 55 | G4cout << "ELYSE::EventAction: cannot get Analysis Tree" << G4endl; |
---|
| 56 | return; |
---|
| 57 | } |
---|
| 58 | //see AIDA::IManagedObject a comment on dynamic_cast |
---|
| 59 | AIDA::IManagedObject* obj; //generic Object managed by a Tree |
---|
| 60 | |
---|
| 61 | obj = usrTree->find("Event"); |
---|
| 62 | if(!obj) { |
---|
| 63 | G4cout << "EventAction: WARNING: no tuple Event" << G4endl; |
---|
| 64 | usrTree->ls(); |
---|
| 65 | //exit(0); |
---|
| 66 | } else { |
---|
| 67 | eventTuple = CAST(obj,AIDA::ITuple); |
---|
| 68 | if (!eventTuple) { |
---|
| 69 | G4cout << "EventAction: FATAL: eventTuple not a Tuple" << G4endl; |
---|
| 70 | usrTree->ls(); |
---|
| 71 | exit(0); |
---|
| 72 | } |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | }//Ctor |
---|
| 76 | |
---|
| 77 | //------------------------------------------------------------------------------------------- |
---|
| 78 | |
---|
| 79 | ELYSE::EventAction::~EventAction(){ |
---|
| 80 | }//Dtor |
---|
| 81 | |
---|
| 82 | //------------------------------------------------------------------------------------------- |
---|
| 83 | |
---|
| 84 | void ELYSE::EventAction::BeginOfEventAction(const G4Event* evt){ |
---|
| 85 | |
---|
| 86 | G4cout << " (JEC) EventAction::Begin EventAction" << G4endl; |
---|
| 87 | if (!evt) evt->Print(); |
---|
| 88 | |
---|
| 89 | }//BeginOfEventAction |
---|
| 90 | |
---|
| 91 | //------------------------------------------------------------------------------------------- |
---|
| 92 | |
---|
| 93 | void ELYSE::EventAction::EndOfEventAction(const G4Event* evt) { |
---|
| 94 | G4cout << " (JEC) EventAction::End EventAction" << G4endl; |
---|
| 95 | |
---|
| 96 | //see Analysis.cxx to get the Tuple variables |
---|
| 97 | |
---|
| 98 | // -------------------- |
---|
| 99 | // Get Event Information |
---|
| 100 | // -------------------- |
---|
| 101 | //JEC FIXME: difference between event_id & vecRecNumber |
---|
| 102 | G4int event_id = evt->GetEventID(); |
---|
| 103 | |
---|
| 104 | G4ThreeVector vtx = generatorAction->GetVtx(); |
---|
| 105 | G4int vtxvol = EventFindStartingVolume(vtx); |
---|
| 106 | |
---|
| 107 | G4cout << ">>>>>>>>>>>>> New Event "; |
---|
| 108 | G4cout << " evt Id " << event_id |
---|
| 109 | << " start in volume " << vtxvol |
---|
| 110 | << G4endl; |
---|
| 111 | G4cout << "Vertex (" |
---|
| 112 | << vtx.x()/cm << " , " |
---|
| 113 | << vtx.y()/cm << " , " |
---|
| 114 | << vtx.z()/cm << " , " |
---|
| 115 | << ")" |
---|
| 116 | << G4endl; |
---|
| 117 | |
---|
| 118 | //JEC FIXME: introduce enumeration for the column shared by Analysis/EventAction & namespace protected |
---|
| 119 | |
---|
| 120 | if(!eventTuple) return; |
---|
| 121 | |
---|
| 122 | eventTuple->fill(0, event_id); //eventId |
---|
| 123 | eventTuple->fill(1, vtxvol); //vtxVol |
---|
| 124 | |
---|
| 125 | AIDA::ITuple* vtxPos = eventTuple->getTuple( 2 ); |
---|
| 126 | vtxPos->fill(0, vtx.x()/cm); //vtxPos |
---|
| 127 | vtxPos->fill(1, vtx.y()/cm); |
---|
| 128 | vtxPos->fill(2, vtx.z()/cm); |
---|
| 129 | vtxPos->addRow(); |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | // -------------------- |
---|
| 133 | // Get Track Information |
---|
| 134 | // -------------------- |
---|
| 135 | |
---|
| 136 | //variables used to fill the tuple (track part) |
---|
| 137 | G4int pId, parent; |
---|
| 138 | G4float timeStart; |
---|
| 139 | G4double dx, dy, dz; |
---|
| 140 | G4double ETot, px, py, pz; |
---|
| 141 | G4int startVol, stopVol; |
---|
| 142 | |
---|
| 143 | //---------------- |
---|
| 144 | // The beam: 1 particle for the moment |
---|
| 145 | //---------------- |
---|
| 146 | AIDA::ITuple* track = eventTuple->getTuple( 4 ); |
---|
| 147 | |
---|
| 148 | G4int beampdg = generatorAction->GetBeamPDG(); |
---|
| 149 | pId = beampdg; //pId |
---|
| 150 | track->fill(0, pId); |
---|
| 151 | |
---|
| 152 | parent = 0; //parent (none) |
---|
| 153 | track->fill(1, parent); |
---|
| 154 | |
---|
| 155 | timeStart = 0; //creation time |
---|
| 156 | track->fill(2, timeStart); |
---|
| 157 | |
---|
| 158 | G4ThreeVector beamdir = generatorAction->GetBeamDir(); //direction |
---|
| 159 | AIDA::ITuple* direction = track->getTuple( 3 ); |
---|
| 160 | dx = beamdir.x(); |
---|
| 161 | dy = beamdir.y(); |
---|
| 162 | dz = beamdir.z(); |
---|
| 163 | |
---|
| 164 | direction->fill(0, dx); |
---|
| 165 | direction->fill(1, dy); |
---|
| 166 | direction->fill(2, dz); |
---|
| 167 | direction->addRow(); |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | G4double beamenergy = generatorAction->GetBeamEnergy(); |
---|
| 171 | G4double mass = G4ParticleTable::GetParticleTable()->FindParticle(beampdg)->GetPDGMass(); |
---|
| 172 | |
---|
| 173 | track->fill(4, mass); |
---|
| 174 | |
---|
| 175 | ETot = beamenergy; |
---|
| 176 | G4double pTot = ETot; |
---|
| 177 | if ( (beampdg != 22) && //gamma & Optical Photon (To be verified) |
---|
| 178 | (std::abs(beampdg) != 12) && //(anti)nu-e |
---|
| 179 | (std::abs(beampdg) != 14) && //(anti)nu-mu |
---|
| 180 | (std::abs(beampdg) != 16) ) { //(anti)nu-tau |
---|
| 181 | pTot = ETot*ETot - mass*mass; |
---|
| 182 | if (pTot < 1e-6){ |
---|
| 183 | pTot = 0.; |
---|
| 184 | } else { |
---|
| 185 | pTot = std::sqrt(pTot); |
---|
| 186 | } |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | track->fill(5, pTot); |
---|
| 190 | track->fill(6, ETot); |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | AIDA::ITuple* momentum = track->getTuple( 7 ); //momentum |
---|
| 194 | px = pTot * dx; |
---|
| 195 | py = pTot * dy; |
---|
| 196 | pz = pTot * dz; |
---|
| 197 | momentum->fill(0, px); |
---|
| 198 | momentum->fill(1, py); |
---|
| 199 | momentum->fill(2, pz); |
---|
| 200 | momentum->addRow(); |
---|
| 201 | |
---|
| 202 | AIDA::ITuple* startPos = track->getTuple( 8 ); //start position |
---|
| 203 | startPos->fill(0, vtx.x()/cm); |
---|
| 204 | startPos->fill(1, vtx.y()/cm); |
---|
| 205 | startPos->fill(2, vtx.z()/cm); |
---|
| 206 | startPos->addRow(); |
---|
| 207 | |
---|
| 208 | AIDA::ITuple* stopPos = track->getTuple( 9 ); //stop position |
---|
| 209 | stopPos->fill(0, vtx.x()/cm); |
---|
| 210 | stopPos->fill(1, vtx.y()/cm); |
---|
| 211 | stopPos->fill(2, vtx.z()/cm); |
---|
| 212 | stopPos->addRow(); |
---|
| 213 | |
---|
| 214 | startVol = stopVol = -1; |
---|
| 215 | track->fill(10, startVol); //Starting Volume |
---|
| 216 | track->fill(11, stopVol); //Stopping Volume |
---|
| 217 | |
---|
| 218 | //add the Beam track to tuple |
---|
| 219 | track->addRow(); |
---|
| 220 | |
---|
| 221 | G4cout << "----> Tk{Beam}: " |
---|
| 222 | << " pId " << pId |
---|
| 223 | << " parent " << parent |
---|
| 224 | << " creation time " << timeStart |
---|
| 225 | << " Volumes " << startVol << " " << stopVol << "\n" |
---|
| 226 | << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n" |
---|
| 227 | << " ETot " << ETot |
---|
| 228 | << " Mass " << mass |
---|
| 229 | << " pTot " << pTot |
---|
| 230 | << " px,py,pz " << px << " " << py << " " << pz << "\n" |
---|
| 231 | << G4endl; |
---|
| 232 | |
---|
| 233 | // -------------------------- |
---|
| 234 | // Loop over Trajectories |
---|
| 235 | // -------------------------- |
---|
| 236 | |
---|
| 237 | G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer(); |
---|
| 238 | |
---|
| 239 | G4int n_trajectories = 0; |
---|
| 240 | if (trajectoryContainer) n_trajectories = trajectoryContainer->entries(); |
---|
| 241 | |
---|
| 242 | G4int trjIndex; //this is the trajectory Index |
---|
| 243 | G4int parentIndex; //this the parent trajectory index |
---|
| 244 | |
---|
| 245 | G4ThreeVector mom; //px,py,pz |
---|
| 246 | G4double mommag; //p = sqrt(px*px + ...) |
---|
| 247 | G4ThreeVector stop; //stopping point |
---|
| 248 | G4ThreeVector start; //starting point |
---|
| 249 | G4String stopVolumeName; //Volume name of the stopping point for retreive |
---|
| 250 | |
---|
| 251 | eventTuple->fill(5, n_trajectories); //nPart |
---|
| 252 | G4cout << "Final # of tracks : " << n_trajectories << G4endl; |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | for (G4int i=0; i < n_trajectories; i++) { |
---|
| 256 | Trajectory* trj = (Trajectory*)((*(evt->GetTrajectoryContainer()))[i]); |
---|
| 257 | |
---|
| 258 | pId = trj->GetPDGEncoding(); //pId |
---|
| 259 | trjIndex = trj->GetTrackID(); |
---|
| 260 | |
---|
| 261 | if ( ! trj->GetSaveFlag() ) continue; //only consider trajectories which are flaged (see TrackingAction, SteppingAction) |
---|
| 262 | |
---|
| 263 | // initial point of the trajectory |
---|
| 264 | G4TrajectoryPoint* aa = (G4TrajectoryPoint*)trj->GetPoint(0) ; |
---|
| 265 | |
---|
| 266 | track->fill(0, pId); //pId |
---|
| 267 | |
---|
| 268 | parentIndex = trj->GetParentID(); |
---|
| 269 | track->fill(1,parentIndex); //parent |
---|
| 270 | |
---|
| 271 | timeStart = trj->GetGlobalTime()/ns; |
---|
| 272 | track->fill(2, timeStart); //creation time |
---|
| 273 | |
---|
| 274 | mom = trj->GetInitialMomentum(); |
---|
| 275 | mommag = mom.mag(); |
---|
| 276 | direction = track->getTuple( 3 ); //direction |
---|
| 277 | dx = mom.x()/mommag; |
---|
| 278 | dy = mom.y()/mommag; |
---|
| 279 | dz = mom.z()/mommag; |
---|
| 280 | direction->fill(0, dx); |
---|
| 281 | direction->fill(1, dy); |
---|
| 282 | direction->fill(2, dz); |
---|
| 283 | direction->addRow(); |
---|
| 284 | |
---|
| 285 | mass = trj->GetParticleDefinition()->GetPDGMass(); |
---|
| 286 | track->fill(4, mass); //mass |
---|
| 287 | |
---|
| 288 | pTot = mommag; |
---|
| 289 | track->fill(5, pTot); //p Total |
---|
| 290 | |
---|
| 291 | ETot = sqrt(mom.mag2() + mass*mass); |
---|
| 292 | track->fill(6, ETot); //Total energy |
---|
| 293 | |
---|
| 294 | momentum = track->getTuple( 7 ); //3-momentum |
---|
| 295 | px = mom.x(); |
---|
| 296 | py = mom.y(); |
---|
| 297 | pz = mom.z(); |
---|
| 298 | momentum->fill(0, px); |
---|
| 299 | momentum->fill(1, py); |
---|
| 300 | momentum->fill(2, pz); |
---|
| 301 | momentum->addRow(); |
---|
| 302 | |
---|
| 303 | start = aa->GetPosition(); |
---|
| 304 | startPos = track->getTuple( 8 ); //start position |
---|
| 305 | startPos->fill(0, start.x()/cm); |
---|
| 306 | startPos->fill(1, start.y()/cm); |
---|
| 307 | startPos->fill(2, start.z()/cm); |
---|
| 308 | startPos->addRow(); |
---|
| 309 | |
---|
| 310 | stop = trj->GetStoppingPoint(); |
---|
| 311 | stopPos = track->getTuple( 9 ); //stop position |
---|
| 312 | stopPos->fill(0, stop.x()/cm); |
---|
| 313 | stopPos->fill(1, stop.y()/cm); |
---|
| 314 | stopPos->fill(2, stop.z()/cm); |
---|
| 315 | stopPos->addRow(); |
---|
| 316 | |
---|
| 317 | startVol = EventFindStartingVolume(start); |
---|
| 318 | track->fill(10, startVol); //Starting Volume |
---|
| 319 | |
---|
| 320 | |
---|
| 321 | stopVolumeName = trj->GetStoppingVolume()->GetName(); |
---|
| 322 | stopVol = EventFindStoppingVolume(stopVolumeName); |
---|
| 323 | track->fill(11, stopVol); //Stopping Volume |
---|
| 324 | |
---|
| 325 | //add the new track to tuple |
---|
| 326 | track->addRow(); |
---|
| 327 | |
---|
| 328 | G4cout << "----> Tk{"<<i<<"}: " |
---|
| 329 | << " pId " << pId |
---|
| 330 | << " parent " << parent |
---|
| 331 | << " creation time " << timeStart |
---|
| 332 | << " Volumes " << startVol << " " << stopVol << "\n" |
---|
| 333 | << " Start Pos (" << start.x()/cm << "," << start.y() << "," << start.z() << ")\n" |
---|
| 334 | << " Stop Pos (" << stop.x()/cm << "," << stop.y() << "," << stop.z() << ")\n" |
---|
| 335 | << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n" |
---|
| 336 | << " m " << mass |
---|
| 337 | << " ETot " << ETot |
---|
| 338 | << " pTot " << pTot |
---|
| 339 | << " px,py,pz " << px << " " << py << " " << pz << "\n" |
---|
| 340 | << G4endl; |
---|
| 341 | |
---|
| 342 | }//end of loop on trajectories |
---|
| 343 | |
---|
| 344 | |
---|
| 345 | |
---|
| 346 | // -------------------- |
---|
| 347 | // Get Hit Collection |
---|
| 348 | // -------------------- |
---|
[294] | 349 | //JEC 22/2/07 not yet available |
---|
[286] | 350 | |
---|
[294] | 351 | // G4SDManager* SDman = G4SDManager::GetSDMpointer(); //JEC FIXME: use data member |
---|
[286] | 352 | |
---|
[294] | 353 | // // Get Hit collection of this event |
---|
| 354 | // G4HCofThisEvent* HCE = evt->GetHCofThisEvent(); |
---|
| 355 | // TrappVolHitsCollection* TrappVolHC = 0; |
---|
[286] | 356 | |
---|
[294] | 357 | // AIDA::ITuple* hit = eventTuple->getTuple(6); //hit |
---|
| 358 | // G4int nHits=0; |
---|
[286] | 359 | |
---|
[294] | 360 | // if (HCE) { |
---|
| 361 | // G4String name = "TrappingVolume"; //see TrappingVolume.cxx Ctor |
---|
| 362 | // G4int collectionID = SDman->GetCollectionID(name); |
---|
| 363 | // TrappVolHC = (TrappVolHitsCollection*)HCE->GetHC(collectionID); |
---|
| 364 | // } |
---|
[286] | 365 | |
---|
[294] | 366 | // if (TrappVolHC) { |
---|
[286] | 367 | |
---|
[294] | 368 | // nHits = TrappVolHC->entries(); |
---|
| 369 | // eventTuple->fill(5, nHits); //nHits |
---|
[286] | 370 | |
---|
[294] | 371 | // TrappVolHit* aHit; |
---|
[286] | 372 | |
---|
[294] | 373 | // for (G4int i=0; i<nHits ;i++) { |
---|
| 374 | // aHit = (*TrappVolHC)[i]; |
---|
[286] | 375 | |
---|
[294] | 376 | // //dump |
---|
| 377 | // aHit->Print(); |
---|
[286] | 378 | |
---|
[294] | 379 | // hit->fill(0,aHit->GetTime()); //time of hit creation |
---|
| 380 | // hit->fill(1,aHit->GetEtot()); //Total energy of the particle creating the hit |
---|
| 381 | // hit->fill(2,aHit->GetCoordinates().x()); //position of the hit |
---|
| 382 | // hit->fill(3,aHit->GetCoordinates().y()); |
---|
| 383 | // hit->fill(4,aHit->GetCoordinates().z()); |
---|
| 384 | // hit->fill(5,aHit->GetTrkId()); //Id of the particle |
---|
| 385 | // hit->fill(6,aHit->GetParentId()); //Id of the parent of the particle |
---|
| 386 | // hit->fill(7,aHit->GetPDGEncoding()); //PDG code of the particle |
---|
| 387 | // hit->fill(8,aHit->GetEdep()); //Energy deposit |
---|
[286] | 388 | |
---|
[294] | 389 | // hit->addRow(); |
---|
| 390 | // }//loop on hits |
---|
| 391 | // }//Hit container |
---|
[286] | 392 | |
---|
| 393 | //Save the Event |
---|
| 394 | eventTuple->addRow(); |
---|
| 395 | |
---|
| 396 | }//EndOfEventAction |
---|
| 397 | |
---|
| 398 | //-------------------------------------------------------------------------------------------------- |
---|
| 399 | |
---|
| 400 | G4int ELYSE::EventAction::EventFindStartingVolume(G4ThreeVector vtx) { |
---|
| 401 | // Get volume of starting point (see GEANT4 FAQ) |
---|
| 402 | |
---|
| 403 | G4Navigator* tmpNavigator = |
---|
| 404 | G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); |
---|
| 405 | |
---|
| 406 | G4VPhysicalVolume* tmpVolume = tmpNavigator->LocateGlobalPointAndSetup(vtx); |
---|
| 407 | G4String vtxVolumeName = tmpVolume->GetName(); |
---|
| 408 | |
---|
| 409 | return FindVolume(vtxVolumeName); |
---|
| 410 | |
---|
| 411 | }//EventFindStartingVolume |
---|
| 412 | |
---|
| 413 | //-------------------------------------------------------------------------------------------------- |
---|
| 414 | |
---|
| 415 | G4int ELYSE::EventAction::EventFindStoppingVolume(G4String stopVolumeName) { |
---|
| 416 | return FindVolume(stopVolumeName); |
---|
| 417 | }//EventFindStoppingVolume |
---|
| 418 | |
---|
| 419 | |
---|
| 420 | //-------------------------------------------------------------------------------------------------- |
---|
| 421 | G4int ELYSE::EventAction::FindVolume(G4String volumeName) { |
---|
| 422 | G4int volumeCode = -1; |
---|
| 423 | |
---|
| 424 | G4cout << "(FindVolume): VolumeName: <" |
---|
| 425 | << volumeName |
---|
| 426 | << ">" |
---|
| 427 | << G4endl; |
---|
| 428 | |
---|
| 429 | if ( volumeName == "expHall" ) { |
---|
| 430 | volumeCode = 0; |
---|
| 431 | } else if ( volumeName == "SensiVol" ) { |
---|
| 432 | volumeCode = 1; |
---|
| 433 | } else if (volumeName == "TyvekSheet") { |
---|
| 434 | volumeCode = 2; |
---|
| 435 | } else if (volumeName == "waterTank") { |
---|
| 436 | volumeCode = 3; |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | return volumeCode; |
---|
| 440 | }//FindVolume |
---|