← Back to branch summary

~c-e-pidcott/maus/devel

« back to all changes in this revision

Viewing changes to src/common_cpp/DetModel/SciFi/SciFiSD.cc

  • Committer: Chris Rogers
  • Date: 2011-11-04 16:46:18 UTC
  • mfrom: (656.1.27 maus-trunk3)
  • Revision ID: chris.rogers@stfc.ac.uk-20111104164618-4gaisprgsgivebix
Changes for release 0.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/DetModel/SciFi/SciFiSD.cc
 
1
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 
2
 *
 
3
 * MAUS is free software: you can redistribute it and/or modify
 
4
 * it under the terms of the GNU General Public License as published by
 
5
 * the Free Software Foundation, either version 3 of the License, or
 
6
 * (at your option) any later version.
 
7
 *
 
8
 * MAUS is distributed in the hope that it will be useful,
 
9
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
 * GNU General Public License for more details.
 
12
 *
 
13
 * You should have received a copy of the GNU General Public License
 
14
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 
15
 *
 
16
 */
 
17
#include "src/common_cpp/DetModel/SciFi/SciFiSD.hh"
 
18
 
 
19
#include <G4TransportationManager.hh>
 
20
#include <G4FieldManager.hh>
 
21
#include <G4Field.hh>
 
22
#include <G4HCofThisEvent.hh>
 
23
#include <G4Step.hh>
 
24
#include <G4ThreeVector.hh>
 
25
#include <G4SDManager.hh>
 
26
#include <G4ios.hh>
 
27
 
 
28
#include <iostream>
 
29
#include <fstream>
 
30
#include <cmath>
 
31
#include "Interface/MICEEvent.hh"
 
32
#include "src/legacy/Config/MiceModule.hh"
 
33
 
 
34
#include "TFile.h"
 
35
#include "TH1F.h"
 
36
#include "TH2F.h"
 
37
#include "TTree.h"
 
38
 
 
39
SciFiSD::SciFiSD(MiceModule* mod) : MAUSSD(mod) {
 
40
}
 
41
 
 
42
 
 
43
G4bool SciFiSD::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) {
 
44
 
 
45
  G4double edep = aStep->GetTotalEnergyDeposit();
 
46
 
 
47
  double pid = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
 
48
 
 
49
  if ( edep == 0. ) return false;
 
50
  // the old chanNo, held for comparison
 
51
  int old_chanNo = legacy_chanNo(aStep);
 
52
 
 
53
  int hit_i = _hits.size();
 
54
  _hits.push_back(Json::Value());
 
55
 
 
56
  G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle();
 
57
  G4int fiberNumber = theTouchable->GetCopyNumber();           // get the fiber copy number
 
58
  G4ThreeVector Pos = aStep->GetPreStepPoint()->GetPosition(); // true MC position
 
59
  G4ThreeVector Mom = aStep->GetTrack()->GetMomentum();        // true momentum
 
60
 
 
61
  Json::Value channel_id;
 
62
  channel_id["type"]           = "Tracker";
 
63
  channel_id["fiber_number"]   = fiberNumber;
 
64
  channel_id["tracker_number"] = _module->propertyInt("Tracker");
 
65
  channel_id["station_number"] = _module->propertyInt("Station");
 
66
  channel_id["plane_number"]   = _module->propertyInt("Plane");
 
67
 
 
68
  _hits[hit_i]["channel_id"] = channel_id;
 
69
  _hits[hit_i]["track_id"]   = aStep->GetTrack()->GetTrackID();
 
70
  _hits[hit_i]["energy"]     = aStep->GetTrack()->GetTotalEnergy();
 
71
  _hits[hit_i]["charge"]     = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
 
72
  _hits[hit_i]["pid"]        = pid;
 
73
  _hits[hit_i]["time"]       = aStep->GetPreStepPoint()->GetGlobalTime();
 
74
  _hits[hit_i]["energy_deposited"] = edep;
 
75
  _hits[hit_i]["momentum"]["x"]    = Mom.x();
 
76
  _hits[hit_i]["momentum"]["y"]    = Mom.y();
 
77
  _hits[hit_i]["momentum"]["z"]    = Mom.z();
 
78
  _hits[hit_i]["position"]["x"]= Pos.x();
 
79
  _hits[hit_i]["position"]["y"]= Pos.y();
 
80
  _hits[hit_i]["position"]["z"]= Pos.z();
 
81
 
 
82
  int chanNo;
 
83
  int numbFibres = 7*2*(_module->propertyDouble("CentralFibre")+0.5);
 
84
  if ( _module->propertyInt("Tracker") == 0 ) {
 
85
    chanNo = floor((numbFibres-fiberNumber)/7);
 
86
  } else {
 
87
    chanNo = floor(fiberNumber/7);
 
88
  }
 
89
 
 
90
  // assert agreement on chanNo with legacy calculation
 
91
  assert(abs(chanNo-old_chanNo) < 2);
 
92
 
 
93
  return true;
 
94
}
 
95
 
 
96
void SciFiSD::EndOfEvent(G4HCofThisEvent* HCE) {
 
97
  // do nothing
 
98
}
 
99
 
 
100
int SciFiSD::legacy_chanNo(G4Step* aStep) {
 
101
  Hep3Vector pos = _module->globalPosition();
 
102
  Hep3Vector perp(-1., 0., 0.);
 
103
  double chanWidth = 1.4945; // effective channel width without overlap
 
104
  perp *= _module->globalRotation();
 
105
  Hep3Vector delta = aStep->GetPreStepPoint()->GetPosition() - pos;
 
106
  double dist = delta.x() * perp.x() + delta.y() * perp.y() + delta.z() * perp.z();
 
107
  double fibre = _module->propertyDouble("CentralFibre") + dist * 2.0 /
 
108
                ( _module->propertyDouble("Pitch") * 7.0 );
 
109
 
 
110
  double centFibre = _module->propertyDouble("CentralFibre");
 
111
  int  firstChan = static_cast<int> (fibre + 0.5);
 
112
  int secondChan(-1);
 
113
  double overlap = chanWidth - (0.1365/2.0);
 
114
  double underlap = (0.1365/2.0);
 
115
  // int nChans(1);
 
116
  nChans = 1;
 
117
  double distInChan = ((firstChan - centFibre)*chanWidth) - dist;
 
118
  // distance in channel greater than section which does not overlap
 
119
  if ((sqrt(distInChan*distInChan) > overlap) || (sqrt(distInChan*distInChan) < underlap)) {
 
120
    nChans = 2;
 
121
    if (distInChan < 0) {
 
122
      secondChan = firstChan - 1;
 
123
    } else {
 
124
      secondChan = firstChan + 1;
 
125
    }
 
126
  }
 
127
  return firstChan;
 
128
}