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