Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file eventgenerator/HepMC/HepMCEx01/src/E 26 /// \file eventgenerator/HepMC/HepMCEx01/src/ExN04StackingAction.cc 27 /// \brief Implementation of the ExN04Stacking 27 /// \brief Implementation of the ExN04StackingAction class 28 // 28 // 29 // << 30 29 31 #include "ExN04StackingAction.hh" 30 #include "ExN04StackingAction.hh" 32 << 31 #include "G4SDManager.hh" 33 #include "ExN04StackingActionMessenger.hh" << 32 #include "G4RunManager.hh" 34 << 35 #include "G4Event.hh" 33 #include "G4Event.hh" 36 #include "G4HCofThisEvent.hh" 34 #include "G4HCofThisEvent.hh" >> 35 #include "G4Track.hh" >> 36 #include "G4TrackStatus.hh" 37 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleTypes.hh" 38 #include "G4ParticleTypes.hh" 39 #include "G4RunManager.hh" << 39 #include "ExN04StackingActionMessenger.hh" 40 #include "G4SDManager.hh" << 41 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 42 #include "G4Track.hh" << 41 #include "G4ios.hh" 43 #include "G4TrackStatus.hh" << 44 42 45 //....oooOO0OOooo........oooOO0OOooo........oo << 46 ExN04StackingAction::ExN04StackingAction() 43 ExN04StackingAction::ExN04StackingAction() 47 : G4UserStackingAction(), fTrkHits(0), fMuon << 44 : trkHits(0), muonHits(0), stage(0) 48 { << 45 { 49 fAngRoI = 30.0 * deg; << 46 angRoI = 30.0*deg; 50 fReqMuon = 2; << 47 reqMuon = 2; 51 fReqIso = 10; << 48 reqIso = 10; 52 fMessenger = new ExN04StackingActionMessenge << 49 theMessenger = new ExN04StackingActionMessenger(this); 53 } 50 } 54 51 55 //....oooOO0OOooo........oooOO0OOooo........oo << 56 ExN04StackingAction::~ExN04StackingAction() 52 ExN04StackingAction::~ExN04StackingAction() 57 { << 53 { delete theMessenger; } 58 delete fMessenger; << 59 } << 60 54 61 //....oooOO0OOooo........oooOO0OOooo........oo << 55 G4ClassificationOfNewTrack 62 G4ClassificationOfNewTrack ExN04StackingAction << 56 ExN04StackingAction::ClassifyNewTrack(const G4Track * aTrack) 63 { 57 { 64 G4ClassificationOfNewTrack classification = 58 G4ClassificationOfNewTrack classification = fWaiting; 65 switch (fStage) { << 59 switch(stage) 66 case 0: // Stage 0 : Primary muons only << 60 { 67 if (aTrack->GetParentID() == 0) { << 61 case 0: // Stage 0 : Primary muons only 68 G4ParticleDefinition* particleType = a << 62 if(aTrack->GetParentID()==0) 69 if ((particleType == G4MuonPlus::MuonP << 63 { 70 || (particleType == G4MuonMinus::M << 64 G4ParticleDefinition * particleType = aTrack->GetDefinition(); 71 { << 65 if((particleType==G4MuonPlus::MuonPlusDefinition()) 72 classification = fUrgent; << 66 ||(particleType==G4MuonMinus::MuonMinusDefinition())) 73 } << 67 { classification = fUrgent; } 74 } << 68 } >> 69 break; >> 70 >> 71 case 1: // Stage 1 : Charged primaries only >> 72 // Suspended tracks will be sent to the waiting stack >> 73 if(aTrack->GetParentID()!=0) { break; } >> 74 if(aTrack->GetTrackStatus()==fSuspend) { break; } >> 75 if(aTrack->GetDefinition()->GetPDGCharge()==0.) { break; } >> 76 classification = fUrgent; >> 77 break; >> 78 >> 79 default: // Stage 2 : Accept all primaries >> 80 // Accept all secondaries in RoI >> 81 // Kill secondaries outside RoI >> 82 if(aTrack->GetParentID()==0) >> 83 { >> 84 classification = fUrgent; 75 break; 85 break; 76 << 86 } 77 case 1: // Stage 1 : Charged primaries on << 87 if((angRoI<0.)||InsideRoI(aTrack,angRoI)) 78 // Suspended tracks wil << 88 { 79 if (aTrack->GetParentID() != 0) { << 80 break; << 81 } << 82 if (aTrack->GetTrackStatus() == fSuspend << 83 break; << 84 } << 85 if (aTrack->GetDefinition()->GetPDGCharg << 86 break; << 87 } << 88 classification = fUrgent; 89 classification = fUrgent; 89 break; 90 break; 90 << 91 } 91 default: // Stage 2 : Accept all primarie << 92 classification = fKill; 92 // Accept all secondar << 93 // Kill secondaries ou << 94 if (aTrack->GetParentID() == 0) { << 95 classification = fUrgent; << 96 break; << 97 } << 98 if ((fAngRoI < 0.) || InsideRoI(aTrack, << 99 classification = fUrgent; << 100 break; << 101 } << 102 classification = fKill; << 103 } 93 } 104 return classification; 94 return classification; 105 } 95 } 106 96 107 //....oooOO0OOooo........oooOO0OOooo........oo << 97 G4bool ExN04StackingAction::InsideRoI(const G4Track * aTrack,G4double ang) 108 G4bool ExN04StackingAction::InsideRoI(const G4 << 109 { 98 { 110 if (!fMuonHits) { << 99 if(!muonHits) 111 fMuonHits = (ExN04MuonHitsCollection*)GetC << 100 { muonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); } 112 } << 101 if(!muonHits) 113 if (!fMuonHits) { << 102 { G4cerr << "muonCollection NOT FOUND" << G4endl; 114 G4cerr << "muonCollection NOT FOUND" << G4 << 103 return true; } 115 return true; << 116 } << 117 104 118 G4int nhits = fMuonHits->entries(); << 105 G4int nhits = muonHits->entries(); 119 106 120 const G4ThreeVector trPos = aTrack->GetPosit 107 const G4ThreeVector trPos = aTrack->GetPosition(); 121 for (G4int i = 0; i < nhits; i++) { << 108 for(G4int i=0;i<nhits;i++) 122 G4ThreeVector muHitPos = (*fMuonHits)[i]-> << 109 { >> 110 G4ThreeVector muHitPos = (*muonHits)[i]->GetPos(); 123 G4double angl = muHitPos.angle(trPos); 111 G4double angl = muHitPos.angle(trPos); 124 if (angl < ang) { << 112 if(angl<ang) { return true; } 125 return true; << 126 } << 127 } 113 } 128 114 129 return false; 115 return false; 130 } 116 } 131 117 132 //....oooOO0OOooo........oooOO0OOooo........oo << 133 G4VHitsCollection* ExN04StackingAction::GetCol 118 G4VHitsCollection* ExN04StackingAction::GetCollection(G4String colName) 134 { 119 { 135 G4SDManager* SDMan = G4SDManager::GetSDMpoin 120 G4SDManager* SDMan = G4SDManager::GetSDMpointer(); 136 G4RunManager* runMan = G4RunManager::GetRunM 121 G4RunManager* runMan = G4RunManager::GetRunManager(); 137 int colID = SDMan->GetCollectionID(colName); 122 int colID = SDMan->GetCollectionID(colName); 138 if (colID >= 0) { << 123 if(colID>=0) >> 124 { 139 const G4Event* currentEvent = runMan->GetC 125 const G4Event* currentEvent = runMan->GetCurrentEvent(); 140 G4HCofThisEvent* HCE = currentEvent->GetHC 126 G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent(); 141 return HCE->GetHC(colID); 127 return HCE->GetHC(colID); 142 } 128 } 143 return 0; 129 return 0; 144 } 130 } 145 131 146 //....oooOO0OOooo........oooOO0OOooo........oo << 147 void ExN04StackingAction::NewStage() 132 void ExN04StackingAction::NewStage() 148 { 133 { 149 fStage++; << 134 stage++; 150 G4int nhits; 135 G4int nhits; 151 if (fStage == 1) { << 136 if(stage==1) 152 // Stage 0->1 : check if at least "fReqMuo << 137 { 153 // otherwise abort current ev << 138 // Stage 0->1 : check if at least "reqMuon" hits on muon chamber 154 if (!fMuonHits) { << 139 // otherwise abort current event 155 fMuonHits = (ExN04MuonHitsCollection*)Ge << 140 if(!muonHits) 156 } << 141 { muonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); } 157 if (!fMuonHits) { << 142 if(!muonHits) 158 G4cerr << "muonCollection NOT FOUND" << << 143 { G4cerr << "muonCollection NOT FOUND" << G4endl; 159 return; << 144 return; } 160 } << 145 nhits = muonHits->entries(); 161 nhits = fMuonHits->entries(); << 146 G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber." 162 G4cout << "Stage 0->1 : " << nhits << " hi << 147 << G4endl; 163 if (nhits < fReqMuon) { << 148 if(nhits<reqMuon) >> 149 { 164 stackManager->clear(); 150 stackManager->clear(); 165 G4cout << "++++++++ event aborted" << G4 151 G4cout << "++++++++ event aborted" << G4endl; 166 return; 152 return; 167 } 153 } 168 stackManager->ReClassify(); 154 stackManager->ReClassify(); 169 return; 155 return; 170 } 156 } 171 157 172 else if (fStage == 2) { << 158 else if(stage==2) 173 // Stage 1->2 : check the isolation of muo << 159 { 174 // at least "fReqIsoMuon" iso << 160 // Stage 1->2 : check the isolation of muon tracks 175 // otherwise abort current ev << 161 // at least "reqIsoMuon" isolated muons 176 // Isolation requires "fReqIs << 162 // otherwise abort current event. 177 // (including own hits) in th << 163 // Isolation requires "reqIso" or less hits 178 // in the tracker layers. << 164 // (including own hits) in the RoI region 179 nhits = fMuonHits->entries(); << 165 // in the tracker layers. 180 if (!fTrkHits) { << 166 nhits = muonHits->entries(); 181 fTrkHits = (ExN04TrackerHitsCollection*) << 167 if(!trkHits) 182 } << 168 { trkHits = (ExN04TrackerHitsCollection*)GetCollection("trackerCollection"); } 183 if (!fTrkHits) { << 169 if(!trkHits) 184 G4cerr << "trackerCollection NOT FOUND" << 170 { G4cerr << "trackerCollection NOT FOUND" << G4endl; 185 return; << 171 return; } 186 } << 172 G4int nTrkhits = trkHits->entries(); 187 G4int nTrkhits = fTrkHits->entries(); << 188 G4int isoMuon = 0; 173 G4int isoMuon = 0; 189 for (G4int j = 0; j < nhits; j++) { << 174 for(G4int j=0;j<nhits;j++) 190 G4ThreeVector hitPos = (*fMuonHits)[j]-> << 175 { >> 176 G4ThreeVector hitPos = (*muonHits)[j]->GetPos(); 191 G4int nhitIn = 0; 177 G4int nhitIn = 0; 192 for (G4int jj = 0; (jj < nTrkhits) && (n << 178 for(G4int jj=0;(jj<nTrkhits)&&(nhitIn<=reqIso);jj++) 193 G4ThreeVector trkhitPos = (*fTrkHits)[ << 179 { 194 if (trkhitPos.angle(hitPos) < fAngRoI) << 180 G4ThreeVector trkhitPos = (*trkHits)[jj]->GetPos(); >> 181 if(trkhitPos.angle(hitPos)<angRoI) nhitIn++; 195 } 182 } 196 if (nhitIn <= fReqIso) isoMuon++; << 183 if(nhitIn<=reqIso) isoMuon++; 197 } 184 } 198 G4cout << "Stage 1->2 : " << isoMuon << " 185 G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl; 199 if (isoMuon < fReqIsoMuon) { << 186 if(isoMuon<reqIsoMuon) >> 187 { 200 stackManager->clear(); 188 stackManager->clear(); 201 G4cout << "++++++++ event aborted" << G4 189 G4cout << "++++++++ event aborted" << G4endl; 202 return; 190 return; 203 } 191 } 204 stackManager->ReClassify(); 192 stackManager->ReClassify(); 205 return; 193 return; 206 } 194 } 207 195 208 else { << 196 else 209 // Other fStage change : just re-classify << 197 { >> 198 // Other stage change : just re-classify 210 stackManager->ReClassify(); 199 stackManager->ReClassify(); 211 } 200 } 212 } 201 } 213 << 202 214 //....oooOO0OOooo........oooOO0OOooo........oo << 215 void ExN04StackingAction::PrepareNewEvent() 203 void ExN04StackingAction::PrepareNewEvent() 216 { << 204 { 217 fStage = 0; << 205 stage = 0; 218 fTrkHits = 0; << 206 trkHits = 0; 219 fMuonHits = 0; << 207 muonHits = 0; 220 } 208 } >> 209 >> 210 221 211