RE05/src/RE05StackingAction.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 /// \file RE05StackingAction.cc
0027 /// \brief Implementation of the RE05StackingAction class
0028
0029 #include "RE05StackingAction.hh"
0030
0031 #include "RE05StackingActionMessenger.hh"
0032
0033 #include "G4Event.hh"
0034 #include "G4HCofThisEvent.hh"
0035 #include "G4ParticleDefinition.hh"
0036 #include "G4ParticleTypes.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4SDManager.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4Track.hh"
0041 #include "G4TrackStatus.hh"
0042 #include "G4ios.hh"
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045
0046 RE05StackingAction::RE05StackingAction()
0047   : G4UserStackingAction(),
0048     fTrkHits(0),
0049     fMuonHits(0),
0050     fMessenger(0),
0051     fStage(0),
0052     fReqMuon(2),
0053     fReqIso(10),
0054     fAngRoI(30.0 * deg)
0055 {
0056   fMessenger = new RE05StackingActionMessenger(this);
0057 }
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 RE05StackingAction::~RE05StackingAction()
0062 {
0063   delete fMessenger;
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 G4ClassificationOfNewTrack RE05StackingAction::ClassifyNewTrack(const G4Track* aTrack)
0069 {
0070   G4ClassificationOfNewTrack classification = fWaiting;
0071   switch (fStage) {
0072     case 0:  // Stage 0 : Primary muons only
0073       if (aTrack->GetParentID() == 0) {
0074         G4ParticleDefinition* particleType = aTrack->GetDefinition();
0075         if ((particleType == G4MuonPlus::MuonPlusDefinition())
0076             || (particleType == G4MuonMinus::MuonMinusDefinition()))
0077         {
0078           classification = fUrgent;
0079         }
0080       }
0081       break;
0082
0083     case 1:  // Stage 1 : Charged primaries only
0084              //           Suspended tracks will be sent to the waiting stack
0085       if (aTrack->GetParentID() != 0) {
0086         break;
0087       }
0088       if (aTrack->GetTrackStatus() == fSuspend) {
0089         break;
0090       }
0091       if (aTrack->GetDefinition()->GetPDGCharge() == 0.) {
0092         break;
0093       }
0094       classification = fUrgent;
0095       break;
0096
0097     default:  // Stage 2 : Accept all primaries
0098               //           Accept all secondaries in RoI
0099               //           Kill secondaries outside RoI
0100       if (aTrack->GetParentID() == 0) {
0101         classification = fUrgent;
0102         break;
0103       }
0104       if ((fAngRoI < 0.) || InsideRoI(aTrack, fAngRoI)) {
0105         classification = fUrgent;
0106         break;
0107       }
0108       classification = fKill;
0109   }
0110   return classification;
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 G4bool RE05StackingAction::InsideRoI(const G4Track* aTrack, G4double ang)
0116 {
0117   if (!fMuonHits) {
0118     fMuonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection");
0119   }
0120   if (!fMuonHits) {
0121     G4cerr << "muonCollection NOT FOUND" << G4endl;
0122     return true;
0123   }
0124
0125   G4int nhits = fMuonHits->entries();
0126
0127   const G4ThreeVector trPos = aTrack->GetPosition();
0128   for (G4int i = 0; i < nhits; i++) {
0129     G4ThreeVector muHitPos = (*fMuonHits)[i]->GetPos();
0130     G4double angl = muHitPos.angle(trPos);
0131     if (angl < ang) {
0132       return true;
0133     }
0134   }
0135
0136   return false;
0137 }
0138
0139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0140
0141 G4VHitsCollection* RE05StackingAction::GetCollection(G4String colName)
0142 {
0143   G4SDManager* SDMan = G4SDManager::GetSDMpointer();
0144   G4RunManager* runMan = G4RunManager::GetRunManager();
0145   int colID = SDMan->GetCollectionID(colName);
0146   if (colID >= 0) {
0147     const G4Event* currentEvent = runMan->GetCurrentEvent();
0148     G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
0149     return HCE->GetHC(colID);
0150   }
0151   return 0;
0152 }
0153
0154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0155
0156 void RE05StackingAction::NewStage()
0157 {
0158   fStage++;
0159   G4int nhits;
0160   if (fStage == 1) {
0161     // Stage 0->1 : check if at least "fReqMuon" hits on muon chamber
0162     //              otherwise abort current event
0163     if (!fMuonHits) {
0164       fMuonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection");
0165     }
0166     if (!fMuonHits) {
0167       G4cerr << "muonCollection NOT FOUND" << G4endl;
0168       return;
0169     }
0170     nhits = fMuonHits->entries();
0171     ////    G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber."
0172     ////         << G4endl;
0173     if (nhits < fReqMuon) {
0174       stackManager->clear();
0175       ////      G4cout << "++++++++ event aborted" << G4endl;
0176       return;
0177     }
0178     stackManager->ReClassify();
0179     return;
0180   }
0181
0182   else if (fStage == 2) {
0183     // Stage 1->2 : check the isolation of muon tracks
0184     //              at least "fReqIsoMuon" isolated muons
0185     //              otherwise abort current event.
0186     //              Isolation requires "fReqIso" or less hits
0187     //              (including own hits) in the RoI region
0188     //              in the tracker layers.
0189     nhits = fMuonHits->entries();
0190     if (!fTrkHits) {
0191       fTrkHits = (RE05TrackerHitsCollection*)GetCollection("trackerCollection");
0192     }
0193     if (!fTrkHits) {
0194       G4cerr << "trackerCollection NOT FOUND" << G4endl;
0195       return;
0196     }
0197     G4int nTrkhits = fTrkHits->entries();
0198     G4int isoMuon = 0;
0199     for (G4int j = 0; j < nhits; j++) {
0200       G4ThreeVector hitPos = (*fMuonHits)[j]->GetPos();
0201       G4int nhitIn = 0;
0202       for (G4int jj = 0; (jj < nTrkhits) && (nhitIn <= fReqIso); jj++) {
0203         G4ThreeVector trkhitPos = (*fTrkHits)[jj]->GetPos();
0204         if (trkhitPos.angle(hitPos) < fAngRoI) nhitIn++;
0205       }
0206       if (nhitIn <= fReqIso) isoMuon++;
0207     }
0208     ////    G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl;
0209     if (isoMuon < fReqIsoMuon) {
0210       stackManager->clear();
0211       ////      G4cout << "++++++++ event aborted" << G4endl;
0212       return;
0213     }
0214     stackManager->ReClassify();
0215     return;
0216   }
0217
0218   else {
0219     // Other fStage change : just re-classify
0220     stackManager->ReClassify();
0221   }
0222 }
0223
0224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0225
0226 void RE05StackingAction::PrepareNewEvent()
0227 {
0228   fStage = 0;
0229   fTrkHits = 0;
0230   fMuonHits = 0;
0231 }
0232
0233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......