TACTICScorer.cc 3.86 KB
Newer Older
1 2
//#define USE_Garfield //only use if compiled with Garfield

Warren's avatar
Warren committed
3 4
#include "TACTICScorer.hh"
#include "G4UnitsTable.hh"
Warren's avatar
Warren committed
5
#include "G4RunManager.hh"
6 7 8 9 10

#ifdef USE_Garfield
#include "GARFDRIFT.h"
#endif

Warren's avatar
Warren committed
11 12
using namespace TACTICScorer;

13
Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int NumberOfSegments, G4int depth, G4double p0, G4double p1, G4double p2, G4double p3) //what do level and depth do?       
Warren's avatar
Warren committed
14 15 16 17 18 19 20 21
:G4VPrimitiveScorer(name, depth),HCID(-1){
  m_ScorerLength = ScorerLength;
  m_NumberOfSegments = NumberOfSegments;
  m_SegmentLength = m_ScorerLength / m_NumberOfSegments;
  m_Level = Level;
  m_Position = G4ThreeVector(-1000,-1000,-1000);
  m_SegmentNumber = -1;
  m_Index = -1;
22 23 24 25
  m_p0 = p0;
  m_p1 = p1;
  m_p2 = p2;
  m_p3 = p3;
Warren's avatar
Warren committed
26 27 28 29 30 31
}

Gas_Scorer::~Gas_Scorer(){}

G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){

32
  G4double* Infos = new G4double[14];
Warren's avatar
Warren committed
33
  m_Position  = aStep->GetPreStepPoint()->GetPosition();
34

Warren's avatar
Warren committed
35
  Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
Warren's avatar
Warren committed
36

Warren's avatar
Warren committed
37 38 39
  Infos[1] = aStep->GetTrack()->GetTrackID();

  Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();;
40
  
Warren's avatar
Warren committed
41 42 43
  Infos[3] = aStep->GetPreStepPoint()->GetGlobalTime();
  Infos[4] = aStep->GetPreStepPoint()->GetKineticEnergy();
  Infos[5] = aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit();
Warren's avatar
Warren committed
44 45

  m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength ) + 1; //Pad number 
Warren's avatar
Warren committed
46 47 48 49 50
  Infos[6] = m_SegmentNumber;
  Infos[7] = m_Position.z();
  Infos[8] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R
  Infos[9] = aStep->GetTrack()->GetVertexPosition()[2];
  Infos[10] = aStep->GetTrack()->GetVertexKineticEnergy();
Warren's avatar
Warren committed
51
  G4ThreeVector p_vec = aStep->GetTrack()->GetVertexMomentumDirection();
Warren's avatar
Warren committed
52 53
  Infos[11] = acos(p_vec[2]/pow(pow(p_vec[0],2)+pow(p_vec[1],2)+pow(p_vec[2],2),0.5))/deg; //angle relative to z axis (theta);   
  Infos[12] = aStep->GetTrack()->GetTrackLength();
54
  Infos[13] = m_p0 + m_p1*Infos[8] + m_p2*Infos[8]*Infos[8] + m_p3*Infos[8]*Infos[8]*Infos[8];
Warren's avatar
Warren committed
55
  
Warren's avatar
Warren committed
56 57
  m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
  m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6;
58
  
Warren's avatar
Warren committed
59
  if(isnan(Infos[10])) {
60 61 62 63
    aStep->GetTrack()->SetTrackStatus(fStopAndKill);
    return 0;
  }
            
64 65 66
#ifdef USE_Garfield
  G4ThreeVector delta_Position = aStep->GetDeltaPosition();

Warren's avatar
Warren committed
67
  GARFDRIFT(Infos[5]/eV, Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0], Infos[2], Infos[11]);
68 69
#endif
  
Warren's avatar
Warren committed
70 71
  map<G4int, G4double**>::iterator it;
  it= EvtMap->GetMap()->find(m_Index);
72
  if(it!=EvtMap->GetMap()->end()){
Warren's avatar
Warren committed
73
    G4double* dummy = *(it->second);
74
    if(Infos[1]==dummy[1]) Infos[5]+=dummy[5]; //accumulate ionisation energy deposit to get total accross pad
Warren's avatar
Warren committed
75 76 77 78
    delete dummy;
  }
  
  EvtMap->set(m_Index, Infos);
79

Warren's avatar
Warren committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
  return TRUE;

}

void Gas_Scorer::Initialize(G4HCofThisEvent* HCE){
    EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
    if (HCID < 0) {
        HCID = GetCollectionID(0);
    }
    HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
}

void Gas_Scorer::EndOfEvent(G4HCofThisEvent*){}

void Gas_Scorer::clear(){
    std::map<G4int, G4double**>::iterator    MapIterator;
    for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
        delete *(MapIterator->second);
    }

    EvtMap->clear();
}


void Gas_Scorer::DrawAll(){}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......                                                                              

void Gas_Scorer::PrintAll(){
    G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl ;
    G4cout << " PrimitiveScorer " << GetName() << G4endl               ;
    G4cout << " Number of entries " << EvtMap->entries() << G4endl     ;
}