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