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