1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.26 2001/05/30 12:18:13 morsch
19 Fastidious printf commented.
21 Revision 1.25 2001/05/23 11:59:46 morsch
22 Use RemoveAt method instead of delete to remove objects from TClonesArray.
24 Revision 1.24 2001/05/21 08:39:13 morsch
25 Use fStepBack = 1 only in debug mode.
27 Revision 1.23 2001/05/20 10:10:39 morsch
28 - Debug output at the beginning of new event and end of run.
29 - Filter out boundary loops.
31 Revision 1.22 2001/05/11 13:22:40 morsch
32 If run with debug option (from gAlice) geantinos are sent back and volume sequence forward/backward is compared.
33 Can be very verbous in some cases.
35 Revision 1.21 2000/12/15 10:33:59 morsch
36 Invert coordinates to make meaningful zylindrical plots.
38 Revision 1.20 2000/11/30 07:12:48 alibrary
39 Introducing new Rndm and QA classes
41 Revision 1.19 2000/10/26 14:13:05 morsch
42 - Change from coordinates theta, phi to general coordinates Coor1 and Coor2.
43 - Lego generator instance can be passed in constructor.
45 Revision 1.18 2000/10/02 21:28:14 fca
46 Removal of useless dependecies via forward declarations
48 Revision 1.17 2000/07/12 08:56:25 fca
49 Coding convention correction and warning removal
51 Revision 1.16 2000/05/26 08:35:03 fca
52 Move the check on z after z has been retrieved
54 Revision 1.15 2000/05/16 13:10:40 fca
55 New method IsNewTrack and fix for a problem in Father-Daughter relations
57 Revision 1.14 2000/04/27 10:38:21 fca
58 Correct termination of Lego Run and introduce Lego getter in AliRun
60 Revision 1.13 2000/04/26 10:17:31 fca
61 Changes in Lego for G4 compatibility
63 Revision 1.12 2000/04/07 11:12:33 fca
64 G4 compatibility changes
66 Revision 1.11 2000/03/22 13:42:26 fca
67 SetGenerator does not replace an existing generator, ResetGenerator does
69 Revision 1.10 2000/02/23 16:25:22 fca
70 AliVMC and AliGeant3 classes introduced
71 ReadEuclid moved from AliRun to AliModule
73 Revision 1.9 1999/12/03 10:54:01 fca
76 Revision 1.8 1999/10/01 09:54:33 fca
77 Correct logics for Lego StepManager
79 Revision 1.7 1999/09/29 09:24:29 fca
80 Introduction of the Copyright and cvs Log
84 //////////////////////////////////////////////////////////////
85 //////////////////////////////////////////////////////////////
87 // Utility class to evaluate the material budget from
88 // a given radius to the surface of an arbitrary cylinder
89 // along radial directions from the centre:
92 // - Interaction length
95 // Geantinos are shot in the bins in the fNtheta bins in theta
96 // and fNphi bins in phi with specified rectangular limits.
97 // The statistics are accumulated per
98 // fRadMin < r < fRadMax and <0 < z < fZMax
100 // To activate this option, you can do:
101 // Root > gAlice.RunLego();
102 // or Root > .x menu.C then select menu item "RunLego"
103 // Note that when calling gAlice->RunLego, an optional list
104 // of arguments may be specified.
108 <img src="picts/alilego.gif">
112 //////////////////////////////////////////////////////////////
118 #include "AliDebugVolume.h"
120 #include "AliLegoGenerator.h"
121 #include "AliConst.h"
124 #include "../TGeant3/TGeant3.h"
126 #include "TClonesArray.h"
131 //___________________________________________
135 // Default constructor
143 //___________________________________________
144 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
145 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
146 Float_t rmin, Float_t rmax, Float_t zmax)
147 : TNamed("Lego Generator",title)
150 // specify the angular limits and the size of the rectangular box
152 fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
153 nphi, phimin, phimax, rmin, rmax, zmax);
157 fHistRadl = new TH2F("hradl","Radiation length map",
158 ntheta,thetamin,thetamax,nphi,phimin,phimax);
159 fHistAbso = new TH2F("habso","Interaction length map",
160 ntheta,thetamin,thetamax,nphi,phimin,phimax);
161 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
162 ntheta,thetamin,thetamax,nphi,phimin,phimax);
168 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
169 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
170 fDebug = gAlice->GetDebug();
174 AliLego::AliLego(const char *title, AliLegoGenerator* generator)
175 : TNamed("Lego Generator",title)
178 // specify the angular limits and the size of the rectangular box
181 Float_t c1min, c1max, c2min, c2max;
182 Int_t n1 = fGener->NCoor1();
183 Int_t n2 = fGener->NCoor2();
185 fGener->Coor1Range(c1min, c1max);
186 fGener->Coor2Range(c2min, c2max);
188 fHistRadl = new TH2F("hradl","Radiation length map",
189 n2, c2min, c2max, n1, c1min, c1max);
190 fHistAbso = new TH2F("habso","Interaction length map",
191 n2, c2min, c2max, n1, c1min, c1max);
192 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
193 n2, c2min, c2max, n1, c1min, c1max);
200 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
201 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
202 fDebug = gAlice->GetDebug();
206 //___________________________________________
216 if (fVolumesFwd) delete fVolumesFwd;
217 if (fVolumesBwd) delete fVolumesBwd;
220 //___________________________________________
221 void AliLego::BeginEvent()
224 // --- Set to 0 radiation length, absorption length and g/cm2 ---
231 if (fErrorCondition) DumpVolumes();
232 fVolumesFwd->Delete();
233 fVolumesBwd->Delete();
237 if (gAlice->CurrentTrack() == 0) fStepBack = 0;
241 //___________________________________________
242 void AliLego::FinishEvent()
245 // Finish the event and update the histos
248 c1 = fGener->CurCoor1();
249 c2 = fGener->CurCoor2();
250 fHistRadl->Fill(c2,c1,fTotRadl);
251 fHistAbso->Fill(c2,c1,fTotAbso);
252 fHistGcm2->Fill(c2,c1,fTotGcm2);
255 //___________________________________________
256 void AliLego::FinishRun()
259 // Store histograms in current Root file
265 // Delete histograms from memory
266 fHistRadl->Delete(); fHistRadl=0;
267 fHistAbso->Delete(); fHistAbso=0;
268 fHistGcm2->Delete(); fHistGcm2=0;
270 if (fErrorCondition) DumpVolumes();
273 //___________________________________________
274 void AliLego::Copy(AliLego &lego) const
277 // Copy *this onto lego -- not implemented
279 Fatal("Copy","Not implemented!\n");
282 //___________________________________________
283 void AliLego::StepManager() {
284 // called from AliRun::Stepmanager from gustep.
285 // Accumulate the 3 parameters step by step
288 Float_t a,z,dens,radl,absl;
291 static Float_t vect[3], dir[3];
295 id = gMC->CurrentVolID(copy);
296 vol = gMC->VolName(id);
297 Float_t step = gMC->TrackStep();
299 TLorentzVector pos, mom;
300 gMC->TrackPosition(pos);
301 gMC->TrackMomentum(mom);
304 if (gMC->IsTrackEntering()) status = 1;
305 if (gMC->IsTrackExiting()) status = 2;
311 // printf("\n volume %s %d", vol, status);
313 // Normal Forward stepping
316 // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
319 // store volume if different from previous
322 TClonesArray &lvols = *fVolumesFwd;
323 if (fStepsForward > 0) {
324 AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesFwd)[fStepsForward-1];
325 if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) {
327 fVolumesFwd->RemoveAt(fStepsForward);
328 fVolumesFwd->RemoveAt(fStepsForward+1);
332 new(lvols[fStepsForward++])
333 AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
337 // Get current material properties
339 gMC->CurrentMaterial(a,z,dens,radl,absl);
343 // --- See if we have to stop now
344 if (TMath::Abs(pos[2]) > fGener->ZMax() ||
345 pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
346 if (!gMC->IsNewTrack()) {
347 // Not the first step, add past contribution
353 // generate "mirror" particle flying back
355 fStepsBackward = fStepsForward;
357 Float_t pmom[3], orig[3];
358 Float_t polar[3] = {0.,0.,0.};
367 gAlice->SetTrack(1, gAlice->CurrentTrack(),
368 0, pmom, orig, polar, 0., kPNoProcess, ntr);
371 } // not a new track !
373 if (fDebug) fStepBack = 1;
376 } // outside scoring region ?
378 // --- See how long we have to go
384 t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
388 fTotAbso += step/absl;
389 fTotRadl += step/radl;
390 fTotGcm2 += step*dens;
397 // Geometry debugging
398 // Fly back and compare volume sequence
400 TClonesArray &lvols = *fVolumesBwd;
401 if (fStepsBackward < fStepsForward) {
402 AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesBwd)[fStepsBackward];
403 if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) {
405 fVolumesBwd->RemoveAt(fStepsBackward-1);
406 fVolumesBwd->RemoveAt(fStepsBackward-2);
411 // printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
412 if (fStepsBackward < 0) {
418 new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
420 AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesFwd)[fStepsBackward];
421 if (! (tmp->IsEqual(vol, copy)) && (!fErrorCondition))
423 printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
424 fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
431 void AliLego::DumpVolumes()
434 // Dump volume sequence in case of error
436 printf("\n Dumping Volume Sequence:");
437 printf("\n ==============================================");
439 for (Int_t i = fStepsForward-1; i>=0; i--)
441 AliDebugVolume* tmp1 = (AliDebugVolume*) (*fVolumesFwd)[i];
442 AliDebugVolume* tmp2 = (AliDebugVolume*) (*fVolumesBwd)[i];
444 printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
446 tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
447 tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
448 if (tmp2 && i>= fStepsBackward)
449 printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
451 tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
452 tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
454 printf("\n ............................................................................\n");
456 printf("\n ==============================================\n");