Corrections to obey our coding conventions
[u/mrichter/AliRoot.git] / STEER / AliLego.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
acd84897 16/* $Id$ */
4c039060 17
fe4da5cc 18//////////////////////////////////////////////////////////////
19//////////////////////////////////////////////////////////////
20//
21// Utility class to evaluate the material budget from
22// a given radius to the surface of an arbitrary cylinder
23// along radial directions from the centre:
24//
25// - radiation length
26// - Interaction length
27// - g/cm2
28//
29// Geantinos are shot in the bins in the fNtheta bins in theta
30// and fNphi bins in phi with specified rectangular limits.
31// The statistics are accumulated per
32// fRadMin < r < fRadMax and <0 < z < fZMax
33//
34// To activate this option, you can do:
35// Root > gAlice.RunLego();
36// or Root > .x menu.C then select menu item "RunLego"
37// Note that when calling gAlice->RunLego, an optional list
38// of arguments may be specified.
39//
40//Begin_Html
41/*
1439f98e 42<img src="picts/alilego.gif">
fe4da5cc 43*/
44//End_Html
45//
46//////////////////////////////////////////////////////////////
47
e2afb3b6 48#include "TClonesArray.h"
49#include "TH2.h"
fe4da5cc 50#include "TMath.h"
e2afb3b6 51#include "TString.h"
5d12ce38 52#include "TVirtualMC.h"
94de3818 53
e2afb3b6 54#include "AliConst.h"
7eb618fa 55#include "AliDebugVolume.h"
e2afb3b6 56#include "AliLego.h"
8918e700 57#include "AliLegoGenerator.h"
e2afb3b6 58#include "AliRun.h"
5d12ce38 59#include "AliMC.h"
7eb618fa 60
fe4da5cc 61ClassImp(AliLego)
62
e2afb3b6 63//_______________________________________________________________________
64AliLego::AliLego():
65 fGener(0),
66 fTotRadl(0),
67 fTotAbso(0),
68 fTotGcm2(0),
69 fHistRadl(0),
70 fHistAbso(0),
71 fHistGcm2(0),
72 fHistReta(0),
73 fVolumesFwd(0),
74 fVolumesBwd(0),
75 fStepBack(0),
76 fStepsBackward(0),
77 fStepsForward(0),
78 fErrorCondition(0),
79 fDebug(0)
fe4da5cc 80{
8918e700 81 //
82 // Default constructor
83 //
fe4da5cc 84}
85
e2afb3b6 86//_______________________________________________________________________
87AliLego::AliLego(const AliLego &lego):
88 TNamed(lego),
89 fGener(0),
90 fTotRadl(0),
91 fTotAbso(0),
92 fTotGcm2(0),
93 fHistRadl(0),
94 fHistAbso(0),
95 fHistGcm2(0),
96 fHistReta(0),
97 fVolumesFwd(0),
98 fVolumesBwd(0),
99 fStepBack(0),
100 fStepsBackward(0),
101 fStepsForward(0),
102 fErrorCondition(0),
103 fDebug(0)
104{
105 //
106 // Copy constructor
107 //
108 lego.Copy(*this);
109}
110
111
112//_______________________________________________________________________
0b27b8c7 113AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6 114 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
115 Float_t rmin, Float_t rmax, Float_t zmax):
116 TNamed("Lego Generator",title),
117 fGener(0),
118 fTotRadl(0),
119 fTotAbso(0),
120 fTotGcm2(0),
121 fHistRadl(0),
122 fHistAbso(0),
123 fHistGcm2(0),
124 fHistReta(0),
125 fVolumesFwd(0),
126 fVolumesBwd(0),
127 fStepBack(0),
128 fStepsBackward(0),
129 fStepsForward(0),
130 fErrorCondition(0),
131 fDebug(0)
fe4da5cc 132{
8918e700 133 //
134 // specify the angular limits and the size of the rectangular box
135 //
0b27b8c7 136 fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077 137 nphi, phimin, phimax, rmin, rmax, zmax);
b13db077 138 fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7 139 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 140 fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7 141 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 142 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7 143 ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa 144//
7eb618fa 145 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
146 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
19237f5f 147 fDebug = gAlice->GetDebug();
0b27b8c7 148}
149
e2afb3b6 150//_______________________________________________________________________
151AliLego::AliLego(const char *title, AliLegoGenerator* generator):
152 TNamed("Lego Generator",title),
153 fGener(0),
154 fTotRadl(0),
155 fTotAbso(0),
156 fTotGcm2(0),
157 fHistRadl(0),
158 fHistAbso(0),
159 fHistGcm2(0),
160 fHistReta(0),
161 fVolumesFwd(0),
162 fVolumesBwd(0),
163 fStepBack(0),
164 fStepsBackward(0),
165 fStepsForward(0),
166 fErrorCondition(0),
167 fDebug(0)
0b27b8c7 168{
169 //
170 // specify the angular limits and the size of the rectangular box
171 //
e2afb3b6 172 fGener = generator;
173 Float_t c1min, c1max, c2min, c2max;
174 Int_t n1 = fGener->NCoor1();
175 Int_t n2 = fGener->NCoor2();
176
177 fGener->Coor1Range(c1min, c1max);
178 fGener->Coor2Range(c2min, c2max);
179
180 fHistRadl = new TH2F("hradl","Radiation length map",
181 n2, c2min, c2max, n1, c1min, c1max);
182 fHistAbso = new TH2F("habso","Interaction length map",
183 n2, c2min, c2max, n1, c1min, c1max);
184 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
185 n2, c2min, c2max, n1, c1min, c1max);
186 //
187 //
b13db077 188
e2afb3b6 189 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
190 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
191 fDebug = gAlice->GetDebug();
fe4da5cc 192}
193
e2afb3b6 194//_______________________________________________________________________
fe4da5cc 195AliLego::~AliLego()
196{
8918e700 197 //
198 // Destructor
199 //
200 delete fHistRadl;
201 delete fHistAbso;
202 delete fHistGcm2;
8918e700 203 delete fGener;
e2afb3b6 204 delete fVolumesFwd;
205 delete fVolumesBwd;
fe4da5cc 206}
207
e2afb3b6 208//_______________________________________________________________________
dffd31ef 209void AliLego::BeginEvent()
b13db077 210{
8918e700 211 //
212 // --- Set to 0 radiation length, absorption length and g/cm2 ---
213 //
dffd31ef 214 fTotRadl = 0;
215 fTotAbso = 0;
216 fTotGcm2 = 0;
19237f5f 217
a74d2699 218 if (fDebug) {
e2afb3b6 219 if (fErrorCondition) DumpVolumes();
220 fVolumesFwd->Delete();
221 fVolumesBwd->Delete();
222 fStepsForward = 0;
223 fStepsBackward = 0;
224 fErrorCondition = 0;
5d12ce38 225 if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
a74d2699 226 }
dffd31ef 227}
228
e2afb3b6 229//_______________________________________________________________________
dffd31ef 230void AliLego::FinishEvent()
231{
8918e700 232 //
233 // Finish the event and update the histos
234 //
0b27b8c7 235 Double_t c1, c2;
236 c1 = fGener->CurCoor1();
237 c2 = fGener->CurCoor2();
10bd7535 238 fHistRadl->Fill(c2,c1,fTotRadl);
239 fHistAbso->Fill(c2,c1,fTotAbso);
240 fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077 241}
242
e2afb3b6 243//_______________________________________________________________________
dffd31ef 244void AliLego::FinishRun()
245{
8918e700 246 //
247 // Store histograms in current Root file
248 //
dffd31ef 249 fHistRadl->Write();
250 fHistAbso->Write();
251 fHistGcm2->Write();
dffd31ef 252
253 // Delete histograms from memory
254 fHistRadl->Delete(); fHistRadl=0;
255 fHistAbso->Delete(); fHistAbso=0;
256 fHistGcm2->Delete(); fHistGcm2=0;
a74d2699 257 //
258 if (fErrorCondition) DumpVolumes();
dffd31ef 259}
260
e2afb3b6 261//_______________________________________________________________________
262void AliLego::Copy(AliLego&) const
8918e700 263{
264 //
265 // Copy *this onto lego -- not implemented
266 //
267 Fatal("Copy","Not implemented!\n");
268}
dffd31ef 269
e2afb3b6 270//_______________________________________________________________________
271void AliLego::StepManager()
272{
273 //
274 // called from AliRun::Stepmanager from gustep.
275 // Accumulate the 3 parameters step by step
276 //
277 static Float_t t;
278 Float_t a,z,dens,radl,absl;
279 Int_t i, id, copy;
280 const char* vol;
281 static Float_t vect[3], dir[3];
282
283 TString tmp1, tmp2;
284 copy = 1;
285 id = gMC->CurrentVolID(copy);
286 vol = gMC->VolName(id);
287 Float_t step = gMC->TrackStep();
288
289 TLorentzVector pos, mom;
290 gMC->TrackPosition(pos);
291 gMC->TrackMomentum(mom);
292
293 Int_t status = 0;
294 if (gMC->IsTrackEntering()) status = 1;
295 if (gMC->IsTrackExiting()) status = 2;
296
7eb618fa 297 if (! fStepBack) {
e2afb3b6 298 // printf("\n volume %s %d", vol, status);
299 //
300 // Normal Forward stepping
301 //
302 if (fDebug) {
303 // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
304
305 //
306 // store volume if different from previous
307 //
7eb618fa 308
a74d2699 309 TClonesArray &lvols = *fVolumesFwd;
7eb618fa 310 if (fStepsForward > 0) {
e2afb3b6 311 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
312 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
a74d2699 313 fStepsForward -= 2;
5f3dc1c9 314 fVolumesFwd->RemoveAt(fStepsForward);
315 fVolumesFwd->RemoveAt(fStepsForward+1);
e2afb3b6 316 }
7eb618fa 317 }
e2afb3b6 318
a74d2699 319 new(lvols[fStepsForward++])
e2afb3b6 320 AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
7eb618fa 321
e2afb3b6 322 } // Debug
323 //
324 // Get current material properties
325
326 gMC->CurrentMaterial(a,z,dens,radl,absl);
327
328 if (z < 1) return;
329
330 // --- See if we have to stop now
331 if (TMath::Abs(pos[2]) > fGener->ZMax() ||
332 pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
333 if (!gMC->IsNewTrack()) {
334 // Not the first step, add past contribution
335 fTotAbso += t/absl;
336 fTotRadl += t/radl;
337 fTotGcm2 += t*dens;
338 if (fDebug) {
339 //
340 // generate "mirror" particle flying back
341 //
342 fStepsBackward = fStepsForward;
343
344 Float_t pmom[3], orig[3];
345 Float_t polar[3] = {0.,0.,0.};
346 Int_t ntr;
347 pmom[0] = -dir[0];
348 pmom[1] = -dir[1];
349 pmom[2] = -dir[2];
350 orig[0] = vect[0];
351 orig[1] = vect[1];
352 orig[2] = vect[2];
353
5d12ce38 354 gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(),
e2afb3b6 355 0, pmom, orig, polar, 0., kPNoProcess, ntr);
356 } // debug
357
358 } // not a new track !
7eb618fa 359
e2afb3b6 360 if (fDebug) fStepBack = 1;
361 gMC->StopTrack();
362 return;
363 } // outside scoring region ?
364
365 // --- See how long we have to go
366 for(i=0;i<3;++i) {
367 vect[i]=pos[i];
368 dir[i]=mom[i];
369 }
370
371 t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
372
373 if(step) {
374
375 fTotAbso += step/absl;
376 fTotRadl += step/radl;
377 fTotGcm2 += step*dens;
378 }
379
7eb618fa 380 } else {
381 if (fDebug) {
e2afb3b6 382 //
383 // Geometry debugging
384 // Fly back and compare volume sequence
385 //
386 TClonesArray &lvols = *fVolumesBwd;
387 if (fStepsBackward < fStepsForward) {
388 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
389 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
390 fStepsBackward += 2;
391 fVolumesBwd->RemoveAt(fStepsBackward-1);
392 fVolumesBwd->RemoveAt(fStepsBackward-2);
7eb618fa 393 }
e2afb3b6 394 }
395
396 fStepsBackward--;
397 // printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
398 if (fStepsBackward < 0) {
7eb618fa 399 gMC->StopTrack();
a74d2699 400 fStepBack = 0;
7eb618fa 401 return;
e2afb3b6 402 }
403
404 new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
405
406 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
407 if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
408 {
409 printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
410 fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
411 fErrorCondition = 1;
412 }
7eb618fa 413 } // Debug
414 } // bwd/fwd
b13db077 415}
416
e2afb3b6 417//_______________________________________________________________________
a74d2699 418void AliLego::DumpVolumes()
419{
e2afb3b6 420 //
421 // Dump volume sequence in case of error
422 //
423 printf("\n Dumping Volume Sequence:");
424 printf("\n ==============================================");
425
426 for (Int_t i = fStepsForward-1; i>=0; i--)
a74d2699 427 {
e2afb3b6 428 AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
429 AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
430 if (tmp1)
431 printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
432 , i,
433 tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
434 tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
435 if (tmp2 && i>= fStepsBackward)
436 printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
437 , i,
438 tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
439 tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
440
441 printf("\n ............................................................................\n");
a74d2699 442 }
e2afb3b6 443 printf("\n ==============================================\n");
a74d2699 444}