Removing AliMCProcess and AliMC
[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
16/*
17$Log$
98490ea9 18Revision 1.30 2002/10/29 14:26:49 hristov
19Code clean-up (F.Carminati)
20
e2afb3b6 21Revision 1.29 2002/10/14 14:57:32 hristov
22Merging the VirtualMC branch to the main development branch (HEAD)
23
b9d0a01d 24Revision 1.28.6.2 2002/07/24 10:08:13 alibrary
25Updating VirtualMC
26
27Revision 1.28.6.1 2002/06/10 14:43:06 hristov
28Merged with v3-08-02
29
30Revision 1.28 2001/10/21 18:38:44 hristov
31Several pointers were set to zero in the default constructors to avoid memory management problems
32
2685bf00 33Revision 1.27 2001/07/20 09:32:18 morsch
34Protection against uncomplete backward stepping in dumping added.
35
899fde02 36Revision 1.26 2001/05/30 12:18:13 morsch
37Fastidious printf commented.
38
09a536f9 39Revision 1.25 2001/05/23 11:59:46 morsch
40Use RemoveAt method instead of delete to remove objects from TClonesArray.
41
5f3dc1c9 42Revision 1.24 2001/05/21 08:39:13 morsch
43Use fStepBack = 1 only in debug mode.
44
19237f5f 45Revision 1.23 2001/05/20 10:10:39 morsch
46- Debug output at the beginning of new event and end of run.
47- Filter out boundary loops.
48
a74d2699 49Revision 1.22 2001/05/11 13:22:40 morsch
50If run with debug option (from gAlice) geantinos are sent back and volume sequence forward/backward is compared.
51Can be very verbous in some cases.
52
7eb618fa 53Revision 1.21 2000/12/15 10:33:59 morsch
54Invert coordinates to make meaningful zylindrical plots.
55
10bd7535 56Revision 1.20 2000/11/30 07:12:48 alibrary
57Introducing new Rndm and QA classes
58
65fb704d 59Revision 1.19 2000/10/26 14:13:05 morsch
60- Change from coordinates theta, phi to general coordinates Coor1 and Coor2.
61- Lego generator instance can be passed in constructor.
62
0b27b8c7 63Revision 1.18 2000/10/02 21:28:14 fca
64Removal of useless dependecies via forward declarations
65
94de3818 66Revision 1.17 2000/07/12 08:56:25 fca
67Coding convention correction and warning removal
68
8918e700 69Revision 1.16 2000/05/26 08:35:03 fca
70Move the check on z after z has been retrieved
71
119a6af6 72Revision 1.15 2000/05/16 13:10:40 fca
73New method IsNewTrack and fix for a problem in Father-Daughter relations
74
a01a8b12 75Revision 1.14 2000/04/27 10:38:21 fca
76Correct termination of Lego Run and introduce Lego getter in AliRun
77
838edcaf 78Revision 1.13 2000/04/26 10:17:31 fca
79Changes in Lego for G4 compatibility
80
dffd31ef 81Revision 1.12 2000/04/07 11:12:33 fca
82G4 compatibility changes
83
875c717b 84Revision 1.11 2000/03/22 13:42:26 fca
85SetGenerator does not replace an existing generator, ResetGenerator does
86
ee1dd322 87Revision 1.10 2000/02/23 16:25:22 fca
88AliVMC and AliGeant3 classes introduced
89ReadEuclid moved from AliRun to AliModule
90
b13db077 91Revision 1.9 1999/12/03 10:54:01 fca
92Fix lego summary
93
00719c1b 94Revision 1.8 1999/10/01 09:54:33 fca
95Correct logics for Lego StepManager
96
f059c84a 97Revision 1.7 1999/09/29 09:24:29 fca
98Introduction of the Copyright and cvs Log
99
4c039060 100*/
101
fe4da5cc 102//////////////////////////////////////////////////////////////
103//////////////////////////////////////////////////////////////
104//
105// Utility class to evaluate the material budget from
106// a given radius to the surface of an arbitrary cylinder
107// along radial directions from the centre:
108//
109// - radiation length
110// - Interaction length
111// - g/cm2
112//
113// Geantinos are shot in the bins in the fNtheta bins in theta
114// and fNphi bins in phi with specified rectangular limits.
115// The statistics are accumulated per
116// fRadMin < r < fRadMax and <0 < z < fZMax
117//
118// To activate this option, you can do:
119// Root > gAlice.RunLego();
120// or Root > .x menu.C then select menu item "RunLego"
121// Note that when calling gAlice->RunLego, an optional list
122// of arguments may be specified.
123//
124//Begin_Html
125/*
1439f98e 126<img src="picts/alilego.gif">
fe4da5cc 127*/
128//End_Html
129//
130//////////////////////////////////////////////////////////////
131
e2afb3b6 132#include "TClonesArray.h"
133#include "TH2.h"
fe4da5cc 134#include "TMath.h"
e2afb3b6 135#include "TString.h"
94de3818 136
e2afb3b6 137#include "AliConst.h"
7eb618fa 138#include "AliDebugVolume.h"
e2afb3b6 139#include "AliLego.h"
8918e700 140#include "AliLegoGenerator.h"
e2afb3b6 141#include "AliRun.h"
7eb618fa 142
fe4da5cc 143ClassImp(AliLego)
144
e2afb3b6 145//_______________________________________________________________________
146AliLego::AliLego():
147 fGener(0),
148 fTotRadl(0),
149 fTotAbso(0),
150 fTotGcm2(0),
151 fHistRadl(0),
152 fHistAbso(0),
153 fHistGcm2(0),
154 fHistReta(0),
155 fVolumesFwd(0),
156 fVolumesBwd(0),
157 fStepBack(0),
158 fStepsBackward(0),
159 fStepsForward(0),
160 fErrorCondition(0),
161 fDebug(0)
fe4da5cc 162{
8918e700 163 //
164 // Default constructor
165 //
fe4da5cc 166}
167
e2afb3b6 168//_______________________________________________________________________
169AliLego::AliLego(const AliLego &lego):
170 TNamed(lego),
171 fGener(0),
172 fTotRadl(0),
173 fTotAbso(0),
174 fTotGcm2(0),
175 fHistRadl(0),
176 fHistAbso(0),
177 fHistGcm2(0),
178 fHistReta(0),
179 fVolumesFwd(0),
180 fVolumesBwd(0),
181 fStepBack(0),
182 fStepsBackward(0),
183 fStepsForward(0),
184 fErrorCondition(0),
185 fDebug(0)
186{
187 //
188 // Copy constructor
189 //
190 lego.Copy(*this);
191}
192
193
194//_______________________________________________________________________
0b27b8c7 195AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6 196 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
197 Float_t rmin, Float_t rmax, Float_t zmax):
198 TNamed("Lego Generator",title),
199 fGener(0),
200 fTotRadl(0),
201 fTotAbso(0),
202 fTotGcm2(0),
203 fHistRadl(0),
204 fHistAbso(0),
205 fHistGcm2(0),
206 fHistReta(0),
207 fVolumesFwd(0),
208 fVolumesBwd(0),
209 fStepBack(0),
210 fStepsBackward(0),
211 fStepsForward(0),
212 fErrorCondition(0),
213 fDebug(0)
fe4da5cc 214{
8918e700 215 //
216 // specify the angular limits and the size of the rectangular box
217 //
0b27b8c7 218 fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077 219 nphi, phimin, phimax, rmin, rmax, zmax);
b13db077 220 fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7 221 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 222 fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7 223 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 224 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7 225 ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa 226//
7eb618fa 227 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
228 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
19237f5f 229 fDebug = gAlice->GetDebug();
0b27b8c7 230}
231
e2afb3b6 232//_______________________________________________________________________
233AliLego::AliLego(const char *title, AliLegoGenerator* generator):
234 TNamed("Lego Generator",title),
235 fGener(0),
236 fTotRadl(0),
237 fTotAbso(0),
238 fTotGcm2(0),
239 fHistRadl(0),
240 fHistAbso(0),
241 fHistGcm2(0),
242 fHistReta(0),
243 fVolumesFwd(0),
244 fVolumesBwd(0),
245 fStepBack(0),
246 fStepsBackward(0),
247 fStepsForward(0),
248 fErrorCondition(0),
249 fDebug(0)
0b27b8c7 250{
251 //
252 // specify the angular limits and the size of the rectangular box
253 //
e2afb3b6 254 fGener = generator;
255 Float_t c1min, c1max, c2min, c2max;
256 Int_t n1 = fGener->NCoor1();
257 Int_t n2 = fGener->NCoor2();
258
259 fGener->Coor1Range(c1min, c1max);
260 fGener->Coor2Range(c2min, c2max);
261
262 fHistRadl = new TH2F("hradl","Radiation length map",
263 n2, c2min, c2max, n1, c1min, c1max);
264 fHistAbso = new TH2F("habso","Interaction length map",
265 n2, c2min, c2max, n1, c1min, c1max);
266 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
267 n2, c2min, c2max, n1, c1min, c1max);
268 //
269 //
b13db077 270
e2afb3b6 271 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
272 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
273 fDebug = gAlice->GetDebug();
fe4da5cc 274}
275
e2afb3b6 276//_______________________________________________________________________
fe4da5cc 277AliLego::~AliLego()
278{
8918e700 279 //
280 // Destructor
281 //
282 delete fHistRadl;
283 delete fHistAbso;
284 delete fHistGcm2;
8918e700 285 delete fGener;
e2afb3b6 286 delete fVolumesFwd;
287 delete fVolumesBwd;
fe4da5cc 288}
289
e2afb3b6 290//_______________________________________________________________________
dffd31ef 291void AliLego::BeginEvent()
b13db077 292{
8918e700 293 //
294 // --- Set to 0 radiation length, absorption length and g/cm2 ---
295 //
dffd31ef 296 fTotRadl = 0;
297 fTotAbso = 0;
298 fTotGcm2 = 0;
19237f5f 299
a74d2699 300 if (fDebug) {
e2afb3b6 301 if (fErrorCondition) DumpVolumes();
302 fVolumesFwd->Delete();
303 fVolumesBwd->Delete();
304 fStepsForward = 0;
305 fStepsBackward = 0;
306 fErrorCondition = 0;
307 if (gAlice->CurrentTrack() == 0) fStepBack = 0;
a74d2699 308 }
dffd31ef 309}
310
e2afb3b6 311//_______________________________________________________________________
dffd31ef 312void AliLego::FinishEvent()
313{
8918e700 314 //
315 // Finish the event and update the histos
316 //
0b27b8c7 317 Double_t c1, c2;
318 c1 = fGener->CurCoor1();
319 c2 = fGener->CurCoor2();
10bd7535 320 fHistRadl->Fill(c2,c1,fTotRadl);
321 fHistAbso->Fill(c2,c1,fTotAbso);
322 fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077 323}
324
e2afb3b6 325//_______________________________________________________________________
dffd31ef 326void AliLego::FinishRun()
327{
8918e700 328 //
329 // Store histograms in current Root file
330 //
dffd31ef 331 fHistRadl->Write();
332 fHistAbso->Write();
333 fHistGcm2->Write();
dffd31ef 334
335 // Delete histograms from memory
336 fHistRadl->Delete(); fHistRadl=0;
337 fHistAbso->Delete(); fHistAbso=0;
338 fHistGcm2->Delete(); fHistGcm2=0;
a74d2699 339 //
340 if (fErrorCondition) DumpVolumes();
dffd31ef 341}
342
e2afb3b6 343//_______________________________________________________________________
344void AliLego::Copy(AliLego&) const
8918e700 345{
346 //
347 // Copy *this onto lego -- not implemented
348 //
349 Fatal("Copy","Not implemented!\n");
350}
dffd31ef 351
e2afb3b6 352//_______________________________________________________________________
353void AliLego::StepManager()
354{
355 //
356 // called from AliRun::Stepmanager from gustep.
357 // Accumulate the 3 parameters step by step
358 //
359 static Float_t t;
360 Float_t a,z,dens,radl,absl;
361 Int_t i, id, copy;
362 const char* vol;
363 static Float_t vect[3], dir[3];
364
365 TString tmp1, tmp2;
366 copy = 1;
367 id = gMC->CurrentVolID(copy);
368 vol = gMC->VolName(id);
369 Float_t step = gMC->TrackStep();
370
371 TLorentzVector pos, mom;
372 gMC->TrackPosition(pos);
373 gMC->TrackMomentum(mom);
374
375 Int_t status = 0;
376 if (gMC->IsTrackEntering()) status = 1;
377 if (gMC->IsTrackExiting()) status = 2;
378
7eb618fa 379 if (! fStepBack) {
e2afb3b6 380 // printf("\n volume %s %d", vol, status);
381 //
382 // Normal Forward stepping
383 //
384 if (fDebug) {
385 // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
386
387 //
388 // store volume if different from previous
389 //
7eb618fa 390
a74d2699 391 TClonesArray &lvols = *fVolumesFwd;
7eb618fa 392 if (fStepsForward > 0) {
e2afb3b6 393 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
394 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
a74d2699 395 fStepsForward -= 2;
5f3dc1c9 396 fVolumesFwd->RemoveAt(fStepsForward);
397 fVolumesFwd->RemoveAt(fStepsForward+1);
e2afb3b6 398 }
7eb618fa 399 }
e2afb3b6 400
a74d2699 401 new(lvols[fStepsForward++])
e2afb3b6 402 AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
7eb618fa 403
e2afb3b6 404 } // Debug
405 //
406 // Get current material properties
407
408 gMC->CurrentMaterial(a,z,dens,radl,absl);
409
410 if (z < 1) return;
411
412 // --- See if we have to stop now
413 if (TMath::Abs(pos[2]) > fGener->ZMax() ||
414 pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
415 if (!gMC->IsNewTrack()) {
416 // Not the first step, add past contribution
417 fTotAbso += t/absl;
418 fTotRadl += t/radl;
419 fTotGcm2 += t*dens;
420 if (fDebug) {
421 //
422 // generate "mirror" particle flying back
423 //
424 fStepsBackward = fStepsForward;
425
426 Float_t pmom[3], orig[3];
427 Float_t polar[3] = {0.,0.,0.};
428 Int_t ntr;
429 pmom[0] = -dir[0];
430 pmom[1] = -dir[1];
431 pmom[2] = -dir[2];
432 orig[0] = vect[0];
433 orig[1] = vect[1];
434 orig[2] = vect[2];
435
436 gAlice->SetTrack(1, gAlice->CurrentTrack(),
437 0, pmom, orig, polar, 0., kPNoProcess, ntr);
438 } // debug
439
440 } // not a new track !
7eb618fa 441
e2afb3b6 442 if (fDebug) fStepBack = 1;
443 gMC->StopTrack();
444 return;
445 } // outside scoring region ?
446
447 // --- See how long we have to go
448 for(i=0;i<3;++i) {
449 vect[i]=pos[i];
450 dir[i]=mom[i];
451 }
452
453 t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
454
455 if(step) {
456
457 fTotAbso += step/absl;
458 fTotRadl += step/radl;
459 fTotGcm2 += step*dens;
460 }
461
7eb618fa 462 } else {
463 if (fDebug) {
e2afb3b6 464 //
465 // Geometry debugging
466 // Fly back and compare volume sequence
467 //
468 TClonesArray &lvols = *fVolumesBwd;
469 if (fStepsBackward < fStepsForward) {
470 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
471 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
472 fStepsBackward += 2;
473 fVolumesBwd->RemoveAt(fStepsBackward-1);
474 fVolumesBwd->RemoveAt(fStepsBackward-2);
7eb618fa 475 }
e2afb3b6 476 }
477
478 fStepsBackward--;
479 // printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
480 if (fStepsBackward < 0) {
7eb618fa 481 gMC->StopTrack();
a74d2699 482 fStepBack = 0;
7eb618fa 483 return;
e2afb3b6 484 }
485
486 new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
487
488 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
489 if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
490 {
491 printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
492 fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
493 fErrorCondition = 1;
494 }
7eb618fa 495 } // Debug
496 } // bwd/fwd
b13db077 497}
498
e2afb3b6 499//_______________________________________________________________________
a74d2699 500void AliLego::DumpVolumes()
501{
e2afb3b6 502 //
503 // Dump volume sequence in case of error
504 //
505 printf("\n Dumping Volume Sequence:");
506 printf("\n ==============================================");
507
508 for (Int_t i = fStepsForward-1; i>=0; i--)
a74d2699 509 {
e2afb3b6 510 AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
511 AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
512 if (tmp1)
513 printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
514 , i,
515 tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
516 tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
517 if (tmp2 && i>= fStepsBackward)
518 printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
519 , i,
520 tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
521 tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
522
523 printf("\n ............................................................................\n");
a74d2699 524 }
e2afb3b6 525 printf("\n ==============================================\n");
a74d2699 526}