3 #include "St_geant_Maker/St_geant_Maker.h"
5 #include "TLorentzVector.h"
6 #include "TGeoManager.h"
7 #include "TGeoManager.h"
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;
20 if ( AgUStep::Instance() ) {
21 (*AgUStep::Instance())();
45 tracks(new TClonesArray(
"Track")),
46 steps(new TClonesArray(
"Step"))
52 Track *Event::AddTrack() {
55 AgUStep::geant3->TrackMomentum( p );
56 AgUStep::geant3->TrackPosition( x );
58 TClonesArray &TRACKS = *tracks;
73 t->mass = AgUStep::geant3->TrackMass();
74 t->charge = AgUStep::geant3->TrackCharge();
79 void Event::Clear(
const Option_t *opts ) {
86 Track::Track() : TObject(), idTruth(-1), eta(0), phi(0), nSteps(0), steps() { Clear(); }
89 Step *Track::AddStep() {
92 AgUStep::geant3->TrackPosition( x );
95 Int_t &n = AgUStep::Instance()->mEvent->nSteps;
96 TClonesArray &STEPS = *AgUStep::Instance()->mEvent->steps;
103 s->idTruth = idTruth;
112 void Track::Clear(
const Option_t *opts ) {
116 px = 0; py = 0; pz = 0;
117 x = 0; y = 0; z = 0; mass = 0; charge = 0;
122 Step::Step() : TObject(),
125 x(0), y(0), z(0), r(0),
139 void Step::Clear(Option_t *opts)
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; }
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);
159 TString Step::volume()
161 TString _volume =
"";
162 if ( gGeoManager )
for ( Int_t i=0;i<15;i++ ) {
164 Int_t
id = vnums[i];
if (0==
id)
break;
165 TGeoVolume *vol = gGeoManager->GetVolume(
id);
166 if ( vol ) _volume = vol->GetName();
174 #define agcstep agcstep_
176 TGiant3 *AgUStep::geant3;
177 Quest_t *AgUStep::cquest;
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;
186 Gcmate_t *AgUStep::cmate;
187 Gccuts_t *AgUStep::ccuts;
188 Gcphys_t *AgUStep::cphys;
189 Gctmed_t *AgUStep::ctmed;
193 AgUStep *AgUStep::sInstance = 0;
195 if ( 0 == sInstance ) sInstance =
new AgUStep();
199 AgUStep::AgUStep() : TNamed(
"AgUStep",
"AgSTAR user stepping routine"),
202 mEvent( new
Event() ),
209 vect0{0.,0.,0., 0.,0.,0., 0.}, oldEvent(-1)
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();
232 void AgUStep::operator()()
235 Double_t x = ctrak -> vect[0];
236 Double_t y = ctrak -> vect[1];
237 Double_t z = ctrak -> vect[2];
239 Double_t _a = cmate->a;
240 Double_t _z = cmate->z;
241 Double_t _dens = cmate->dens;
243 Double_t r = TMath::Sqrt(x*x+y*y);
244 if (r > rmax || TMath::Abs(z) > 200 )
251 idEvent = geant3->CurrentEvent();
254 if ( oldEvent != idEvent )
257 if (mTree && idEvent>1) mTree->Fill();
259 mEvent -> Clear(
"C");
265 mEvent -> idEvent = idEvent;
270 if ( 0 == ctrak->sleng )
279 mTrack = mEvent->AddTrack();
280 mTrack->idTruth = idTruth;
286 aDeStep += ctrak->destep;
287 aStep += ctrak->step;
293 if ( r < rmin )
return;
294 if ( z < zmin || z > zmax )
return;
296 Step *mStep = mTrack -> AddStep();
298 mStep -> isvol = ctmed->isvol;
299 if ( ctmed->isvol ) {
303 mStep->isvol = csets->numbv[0] + 1;
310 mStep -> adEstep = aDeStep;
311 mStep -> nStep = nStep;
312 mStep -> dEstep = ctrak -> destep;
313 mStep -> step = ctrak -> step;
314 mStep -> state = ctrak->inwvol;
318 mStep -> dens = _dens;
326 for ( Int_t i=0;i<cvolu->nlevel;i++ )
330 memcpy( buff, &cvolu->names[i],
sizeof(cvolu->names[i]) );
332 TString volume;
for ( Int_t ii=0;ii<4;ii++ ) volume += buff[ii];
334 path += volume; path +=
"_";
335 path += cvolu->number[i];
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;
357 if ( mTrack->idTruth < mnTruth || mTrack->idTruth > mxTruth )
return;
361 LOG_INFO << Form(
"[AgUStep] idtruth=%i %8.4f %8.4f %8.4f %8.4f %4s %4i %8.4f %8.4f %s",
386 mFile = TFile::Open( filename,
"recreate" );
387 mTree =
new TTree(
"stepping",
"custom stepping tree" );
388 mTree->Branch(
"Event", &mEvent, 32000, 99);
391 void AgUStep::Finish()
void Init(const Char_t *filename="")
Initialize stepping routine. Opens TFile and creates TTree.
TGiant3 * Geant3()
Returns a pointer to the GEANT3 VMC interface.