ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExN04StackingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ExN04StackingAction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 
31 #include "G4SDManager.hh"
32 #include "G4RunManager.hh"
33 #include "G4Event.hh"
34 #include "G4HCofThisEvent.hh"
35 #include "G4Track.hh"
36 #include "G4TrackStatus.hh"
37 #include "G4ParticleDefinition.hh"
38 #include "G4ParticleTypes.hh"
39 #include "ExN04StackingActionMessenger.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "ExN04StackingAction.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46  fTrkHits(0), fMuonHits(0), fStage(0)
47 {
48  fAngRoI = 30.0*deg;
49  fReqMuon = 2;
50  fReqIso = 10;
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 { delete fMessenger; }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 {
62  G4ClassificationOfNewTrack classification = fWaiting;
63  switch(fStage)
64  {
65  case 0: // Stage 0 : Primary muons only
66  if(aTrack->GetParentID()==0)
67  {
68  G4ParticleDefinition * particleType = aTrack->GetDefinition();
69  if((particleType==G4MuonPlus::MuonPlusDefinition())
70  ||(particleType==G4MuonMinus::MuonMinusDefinition()))
71  { classification = fUrgent; }
72  }
73  break;
74 
75  case 1: // Stage 1 : Charged primaries only
76  // Suspended tracks will be sent to the waiting stack
77  if(aTrack->GetParentID()!=0) { break; }
78  if(aTrack->GetTrackStatus()==fSuspend) { break; }
79  if(aTrack->GetDefinition()->GetPDGCharge()==0.) { break; }
80  classification = fUrgent;
81  break;
82 
83  default: // Stage 2 : Accept all primaries
84  // Accept all secondaries in RoI
85  // Kill secondaries outside RoI
86  if(aTrack->GetParentID()==0)
87  {
88  classification = fUrgent;
89  break;
90  }
91  if((fAngRoI<0.)||InsideRoI(aTrack,fAngRoI))
92  {
93  classification = fUrgent;
94  break;
95  }
96  classification = fKill;
97  }
98  return classification;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 {
104  if(!fMuonHits)
105  { fMuonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); }
106  if(!fMuonHits)
107  { G4cerr << "muonCollection NOT FOUND" << G4endl;
108  return true; }
109 
110  G4int nhits = fMuonHits->entries();
111 
112  const G4ThreeVector trPos = aTrack->GetPosition();
113  for(G4int i=0;i<nhits;i++)
114  {
115  G4ThreeVector muHitPos = (*fMuonHits)[i]->GetPos();
116  G4double angl = muHitPos.angle(trPos);
117  if(angl<ang) { return true; }
118  }
119 
120  return false;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 {
128  int colID = SDMan->GetCollectionID(colName);
129  if(colID>=0)
130  {
131  const G4Event* currentEvent = runMan->GetCurrentEvent();
132  G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
133  return HCE->GetHC(colID);
134  }
135  return 0;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 {
141  fStage++;
142  G4int nhits;
143  if(fStage==1)
144  {
145  // Stage 0->1 : check if at least "fReqMuon" hits on muon chamber
146  // otherwise abort current event
147  if(!fMuonHits)
148  { fMuonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); }
149  if(!fMuonHits)
150  { G4cerr << "muonCollection NOT FOUND" << G4endl;
151  return; }
152  nhits = fMuonHits->entries();
153  G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber."
154  << G4endl;
155  if(nhits<fReqMuon)
156  {
157  stackManager->clear();
158  G4cout << "++++++++ event aborted" << G4endl;
159  return;
160  }
162  return;
163  }
164 
165  else if(fStage==2)
166  {
167  // Stage 1->2 : check the isolation of muon tracks
168  // at least "fReqIsoMuon" isolated muons
169  // otherwise abort current event.
170  // Isolation requires "fReqIso" or less hits
171  // (including own hits) in the RoI region
172  // in the tracker layers.
173  nhits = fMuonHits->entries();
174  if(!fTrkHits)
175  { fTrkHits =
176  (ExN04TrackerHitsCollection*)GetCollection("trackerCollection"); }
177  if(!fTrkHits)
178  { G4cerr << "trackerCollection NOT FOUND" << G4endl;
179  return; }
180  G4int nTrkhits = fTrkHits->entries();
181  G4int isoMuon = 0;
182  for(G4int j=0;j<nhits;j++)
183  {
184  G4ThreeVector hitPos = (*fMuonHits)[j]->GetPos();
185  G4int nhitIn = 0;
186  for(G4int jj=0;(jj<nTrkhits)&&(nhitIn<=fReqIso);jj++)
187  {
188  G4ThreeVector trkhitPos = (*fTrkHits)[jj]->GetPos();
189  if(trkhitPos.angle(hitPos)<fAngRoI) nhitIn++;
190  }
191  if(nhitIn<=fReqIso) isoMuon++;
192  }
193  G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl;
194  if(isoMuon<fReqIsoMuon)
195  {
196  stackManager->clear();
197  G4cout << "++++++++ event aborted" << G4endl;
198  return;
199  }
201  return;
202  }
203 
204  else
205  {
206  // Other fStage change : just re-classify
208  }
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 {
214  fStage = 0;
215  fTrkHits = 0;
216  fMuonHits = 0;
217 }