StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AgUStep.cxx
1 #include "AgUStep.h"
2 #include "StMessMgr.h"
3 #include "St_geant_Maker/St_geant_Maker.h"
4 #include "TMath.h"
5 #include "TLorentzVector.h"
6 #include "TGeoManager.h"
7 #include "TGeoManager.h"
8 
9 Float_t AgUStep::rmin = 0.0;
10 Float_t AgUStep::rmax = 200.0;
11 Float_t AgUStep::zmin =-200.0;
12 Float_t AgUStep::zmax = 200.0;
13 Int_t AgUStep::verbose = 0;
14 Int_t AgUStep::mnTruth = 0;
15 Int_t AgUStep::mxTruth = -1;
16 
17 extern "C" {
18 
19  void agustep_() {
20  if ( AgUStep::Instance() ) {
21  (*AgUStep::Instance())();
22  }
23  };
24 
25 
26  struct Agcstep_t {
27 
28  Float_t vect0[7];
29  Float_t vloc0[7];
30  Float_t vloc[7];
31  Float_t xloc[7];
32  Float_t astep;
33  Float_t adestep;
34 
35  } agcstep_;
36 
37 };
38 
39 // .................................................................. Event ...............
40 Event::Event() :
41  TObject(),
42  idEvent(-1),
43  nTracks(0),
44  nSteps(0),
45  tracks(new TClonesArray("Track")),
46  steps(new TClonesArray("Step"))
47 {
48  Clear();
49 }
50 //
51 
52 Track *Event::AddTrack() {
53 
54  TLorentzVector p, x;
55  AgUStep::geant3->TrackMomentum( p );
56  AgUStep::geant3->TrackPosition( x );
57 
58  TClonesArray &TRACKS = *tracks;
59 
60  Track *t = new (TRACKS[nTracks++]) Track();
61 
62  t->eta = p.Eta();
63  t->phi = p.Phi();
64 
65  t->px = p[0];
66  t->py = p[1];
67  t->pz = p[2];
68 
69  t->x = x[0];
70  t->y = x[1];
71  t->z = x[2];
72 
73  t->mass = AgUStep::geant3->TrackMass();
74  t->charge = AgUStep::geant3->TrackCharge();
75 
76  // tracks.Add( t );
77  return t;
78 };
79 void Event::Clear( const Option_t *opts ) {
80  idEvent = -1;
81  nTracks = 0;
82  nSteps = 0;
83  tracks->Clear("C");
84 };
85 // .................................................................. Track ...............
86 Track::Track() : TObject(), idTruth(-1), eta(0), phi(0), nSteps(0), steps() { Clear(); }
87 
88 
89 Step *Track::AddStep() {
90 
91  TLorentzVector x;
92  AgUStep::geant3->TrackPosition( x );
93  // nSteps++;
94  // Step *s = new Step();
95  Int_t &n = AgUStep::Instance()->mEvent->nSteps;
96  TClonesArray &STEPS = *AgUStep::Instance()->mEvent->steps;
97  Step *s = new (STEPS[n++]) Step();
98  s->x = x[0];
99  s->y = x[1];
100  s->z = x[2];
101  s->r = x.Perp();
102  s->idStep = n; // ID of tracking step
103  s->idTruth = idTruth;
104  this->nSteps++; // Why is this not incremented in the track object?
105  steps.Add(s);
106 
107  // LOG_INFO << " New Step at index: " << n << " nSteps=" << nSteps << endm;
108 
109  // nSteps++; // add step
110  return s;
111 }
112 void Track::Clear( const Option_t *opts ) {
113  idTruth = -1;
114  eta = 0;
115  phi = 0;
116  px = 0; py = 0; pz = 0;
117  x = 0; y = 0; z = 0; mass = 0; charge = 0;
118  nSteps=0;
119  steps.Clear("");
120 };
121 // .................................................................. Step ................
122 Step::Step() : TObject(),
123  idStep(-1),
124  idTruth(0),
125  x(0), y(0), z(0), r(0),
126  state(0),
127  dEstep(-1),
128  adEstep(-1),
129  nStep(-1) ,
130  step(-1) ,
131  dens(0) ,
132  A(0),
133  Z(0),
134  isvol(0)
135 {
136  Clear();
137 }
138 
139 void Step::Clear(Option_t *opts)
140 {
141  idStep=-1; idTruth=0; isvol=0;
142  x=0; y=0; z=0; dEstep=-1; adEstep=-1; step=-1; state=0; nStep=-1;
143  for( Int_t i=0;i<15;i++ ) { vnums[i]=0; cnums[i]=0; }
144 }
145 
146 TString Step::path()
147 {
148  TString _path = "";
149  if ( gGeoManager ) for ( Int_t i=0;i<15;i++ ) {
150  Int_t vid = vnums[i]; if (0==vid) break;
151  Int_t cid = cnums[i];
152  TGeoVolume *vol = gGeoManager->GetVolume(vid); if (0==vol) {_path=""; break; }
153  TString name = vol->GetName();
154  _path += Form("/%s_%i",name.Data(),cid);
155  }
156  return _path;
157 }
158 
159 TString Step::volume()
160 {
161  TString _volume = "";
162  if ( gGeoManager ) for ( Int_t i=0;i<15;i++ ) {
163 
164  Int_t id = vnums[i]; if (0==id) break;
165  TGeoVolume *vol = gGeoManager->GetVolume(id);
166  if ( vol ) _volume = vol->GetName();
167  }
168 
169 
170  return _volume;
171 }
172 
173 // .....................................................................................
174 #define agcstep agcstep_
175 
176 TGiant3 *AgUStep::geant3; // G3 VMC
177 Quest_t *AgUStep::cquest; // G3 Commons ...
178 Gclink_t *AgUStep::clink;
179 Gcflag_t *AgUStep::cflag;
180 Gcvolu_t *AgUStep::cvolu;
181 Gcnum_t *AgUStep::cnum;
182 Gcsets_t *AgUStep::csets;
183 Gckine_t *AgUStep::ckine;
184 Gcking_t *AgUStep::cking;
185 Gctrak_t *AgUStep::ctrak;
186 Gcmate_t *AgUStep::cmate;
187 Gccuts_t *AgUStep::ccuts;
188 Gcphys_t *AgUStep::cphys;
189 Gctmed_t *AgUStep::ctmed;
190 
191 Int_t AgUStep::nlev;
192 
193 AgUStep *AgUStep::sInstance = 0;
194 AgUStep *AgUStep::Instance() {
195  if ( 0 == sInstance ) sInstance = new AgUStep();
196  return sInstance;
197 }
198 
199 AgUStep::AgUStep() : TNamed("AgUStep","AgSTAR user stepping routine"),
200  mTree(0),
201  mFile(0),
202  mEvent( new Event() ),
203  mTrack(0),
204  idEvent(-1),
205  idTruth(0),
206  aDeStep(0.),
207  nStep(0),
208  aStep(0.),
209  vect0{0.,0.,0., 0.,0.,0., 0.}, oldEvent(-1)
210 
211 {
212  geant3 = St_geant_Maker::instance()->Geant3();
213  cquest = (Quest_t *) geant3->Quest();
214  clink = (Gclink_t *) geant3->Gclink();
215  cflag = (Gcflag_t *) geant3->Gcflag();
216  cvolu = (Gcvolu_t *) geant3->Gcvolu();
217  cnum = (Gcnum_t *) geant3->Gcnum();
218  csets = (Gcsets_t *) geant3->Gcsets();
219  ckine = (Gckine_t *) geant3->Gckine();
220  cking = (Gcking_t *) geant3->Gcking();
221  ctrak = (Gctrak_t *) geant3->Gctrak();
222  cmate = (Gcmate_t *) geant3->Gcmate();
223  ccuts = (Gccuts_t *) geant3->Gccuts();
224  cphys = (Gcphys_t *) geant3->Gcphys();
225  ctmed = (Gctmed_t *) geant3->Gctmed();
226 
227  oldEvent = -999;
228 
229 };
230 
231 // Take a step through the G3 geometry
232 void AgUStep::operator()()
233 {
234 
235  Double_t x = ctrak -> vect[0];
236  Double_t y = ctrak -> vect[1];
237  Double_t z = ctrak -> vect[2];
238 
239  Double_t _a = cmate->a;
240  Double_t _z = cmate->z;
241  Double_t _dens = cmate->dens;
242 
243  Double_t r = TMath::Sqrt(x*x+y*y);
244  if (r > rmax || TMath::Abs(z) > 200 ) // limited to inner tracking region
245  {
246  ctrak->istop = 2; // stop the track
247  return; // track is exiting region of interest
248  }
249 
250  // Get current event
251  idEvent = geant3->CurrentEvent();
252 
253  // Detect new event and reset sums, clear event, etc...
254  if ( oldEvent != idEvent )
255  {
256 
257  if (mTree && idEvent>1) mTree->Fill(); // last event should be filled on finish
258 
259  mEvent -> Clear("C"); // clear old event
260  aDeStep = 0; // clear sums etc...
261  aStep = 0;
262  nStep = 0;
263  idTruth = 0;
264  // LOG_INFO << "New Event " << idEvent << endm;
265  mEvent -> idEvent = idEvent; // and set the event number
266  oldEvent = idEvent;
267  }
268 
269  // Detect a new track
270  if ( 0 == ctrak->sleng )
271  {
272  aDeStep = 0;
273  aStep = 0;
274  nStep = 0;
275  idTruth++;
276  // LOG_INFO << " New Track " << idTruth << endm;
277 
278  // Add track to this event
279  mTrack = mEvent->AddTrack();
280  mTrack->idTruth = idTruth;
281 
282  }
283 
284 
285 
286  aDeStep += ctrak->destep;
287  aStep += ctrak->step;
288  nStep ++;
289 
290  // Add a new step to the track
291  assert(mTrack);
292 
293  if ( r < rmin ) return; // track has not entered region of interest
294  if ( z < zmin || z > zmax ) return; // outside of region of interest
295 
296  Step *mStep = mTrack -> AddStep();
297 
298  mStep -> isvol = ctmed->isvol;
299  if ( ctmed->isvol ) {
300  // I cannot rely on this, but we will set isvol to the value
301  // of csets->numbv[0]... only correct for sensitive volumes
302  // placed w/in the mother volume ...
303  mStep->isvol = csets->numbv[0] + 1;
304 
305  // LOG_INFO << "Sensitive volume: " << endm;
306  // for ( int i=0;i<20;i++ )
307  // LOG_INFO << csets->numbv[i] << endm;
308  }
309 
310  mStep -> adEstep = aDeStep;
311  mStep -> nStep = nStep;
312  mStep -> dEstep = ctrak -> destep;
313  mStep -> step = ctrak -> step;
314  mStep -> state = ctrak->inwvol;
315 
316  mStep -> A = _a;
317  mStep -> Z = _z;
318  mStep -> dens = _dens;
319 
320  // if (!gGeoManager) return; // step through before complete init?
321 
322  // Print out current path...
323  //LOG_INFO << "N level = " << cvolu->nlevel << endm;
324  TString path = "";
325  TString leaf = "";
326  for ( Int_t i=0;i<cvolu->nlevel;i++ )
327  {
328  path += "/";
329  Char_t buff[16];
330  memcpy( buff, &cvolu->names[i], sizeof(cvolu->names[i]) );
331 
332  TString volume; for ( Int_t ii=0;ii<4;ii++ ) volume += buff[ii];
333 
334  path += volume; path += "_";
335  path += cvolu->number[i];
336 
337  if ( gGeoManager ) {
338  UShort_t volumeNumber = gGeoManager->FindVolumeFast(volume)->GetNumber();
339  UShort_t copyNumber = cvolu->number[i];
340  mStep->vnums[i] = volumeNumber;
341  mStep->cnums[i] = copyNumber;
342  }
343  else {
344  mStep->vnums[i] = 0;
345  mStep->cnums[i] = 0;
346  }
347  leaf = volume;
348 
349  }
350 
352  //
353  // At this point the step is considered ended. We only handle verbose
354  // level printouts below.
355  //
357  if ( mTrack->idTruth < mnTruth || mTrack->idTruth > mxTruth ) return;
358 
359 
360  if ( verbose ) {
361  LOG_INFO << Form("[AgUStep] idtruth=%i %8.4f %8.4f %8.4f %8.4f %4s %4i %8.4f %8.4f %s",
362  mTrack->idTruth,
363  mStep->x,
364  mStep->y,
365  mStep->z,
366  mStep->r,
367  leaf.Data(),
368  mStep->idTruth,
369  aStep,
370  mStep->step,
371  path.Data()
372  )
373  << endm;
374 
375  }
376 
377 }
378 
379 
380 
381 
382 // Initialization
383 void AgUStep::Init( const Char_t *filename )
384 {
385  // LOG_INFO << "Initialize stepper with filename= " << filename << endm;
386  mFile = TFile::Open( filename, "recreate" );
387  mTree = new TTree( "stepping", "custom stepping tree" );
388  mTree->Branch("Event", &mEvent, 32000, 99); // Add the event to the ttree
389 }
390 
391 void AgUStep::Finish()
392 {
393  if (mTree&&mFile) {
394  mTree->Fill();
395 
396  mFile->cd();
397  mTree->Write();
398  // mFile->Write();
399  delete mFile;
400  }
401 }
void Init(const Char_t *filename="")
Initialize stepping routine. Opens TFile and creates TTree.
Definition: AgUStep.cxx:383
TGiant3 * Geant3()
Returns a pointer to the GEANT3 VMC interface.
Definition: AgUStep.h:26
C++ STL includes.
Definition: AgUStep.h:47
Definition: AgUStep.h:71