Replacing Header with Id
[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"
94de3818 52
e2afb3b6 53#include "AliConst.h"
7eb618fa 54#include "AliDebugVolume.h"
e2afb3b6 55#include "AliLego.h"
8918e700 56#include "AliLegoGenerator.h"
e2afb3b6 57#include "AliRun.h"
7eb618fa 58
fe4da5cc 59ClassImp(AliLego)
60
e2afb3b6 61//_______________________________________________________________________
62AliLego::AliLego():
63 fGener(0),
64 fTotRadl(0),
65 fTotAbso(0),
66 fTotGcm2(0),
67 fHistRadl(0),
68 fHistAbso(0),
69 fHistGcm2(0),
70 fHistReta(0),
71 fVolumesFwd(0),
72 fVolumesBwd(0),
73 fStepBack(0),
74 fStepsBackward(0),
75 fStepsForward(0),
76 fErrorCondition(0),
77 fDebug(0)
fe4da5cc 78{
8918e700 79 //
80 // Default constructor
81 //
fe4da5cc 82}
83
e2afb3b6 84//_______________________________________________________________________
85AliLego::AliLego(const AliLego &lego):
86 TNamed(lego),
87 fGener(0),
88 fTotRadl(0),
89 fTotAbso(0),
90 fTotGcm2(0),
91 fHistRadl(0),
92 fHistAbso(0),
93 fHistGcm2(0),
94 fHistReta(0),
95 fVolumesFwd(0),
96 fVolumesBwd(0),
97 fStepBack(0),
98 fStepsBackward(0),
99 fStepsForward(0),
100 fErrorCondition(0),
101 fDebug(0)
102{
103 //
104 // Copy constructor
105 //
106 lego.Copy(*this);
107}
108
109
110//_______________________________________________________________________
0b27b8c7 111AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6 112 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
113 Float_t rmin, Float_t rmax, Float_t zmax):
114 TNamed("Lego Generator",title),
115 fGener(0),
116 fTotRadl(0),
117 fTotAbso(0),
118 fTotGcm2(0),
119 fHistRadl(0),
120 fHistAbso(0),
121 fHistGcm2(0),
122 fHistReta(0),
123 fVolumesFwd(0),
124 fVolumesBwd(0),
125 fStepBack(0),
126 fStepsBackward(0),
127 fStepsForward(0),
128 fErrorCondition(0),
129 fDebug(0)
fe4da5cc 130{
8918e700 131 //
132 // specify the angular limits and the size of the rectangular box
133 //
0b27b8c7 134 fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077 135 nphi, phimin, phimax, rmin, rmax, zmax);
b13db077 136 fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7 137 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 138 fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7 139 ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077 140 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7 141 ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa 142//
7eb618fa 143 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
144 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
19237f5f 145 fDebug = gAlice->GetDebug();
0b27b8c7 146}
147
e2afb3b6 148//_______________________________________________________________________
149AliLego::AliLego(const char *title, AliLegoGenerator* generator):
150 TNamed("Lego Generator",title),
151 fGener(0),
152 fTotRadl(0),
153 fTotAbso(0),
154 fTotGcm2(0),
155 fHistRadl(0),
156 fHistAbso(0),
157 fHistGcm2(0),
158 fHistReta(0),
159 fVolumesFwd(0),
160 fVolumesBwd(0),
161 fStepBack(0),
162 fStepsBackward(0),
163 fStepsForward(0),
164 fErrorCondition(0),
165 fDebug(0)
0b27b8c7 166{
167 //
168 // specify the angular limits and the size of the rectangular box
169 //
e2afb3b6 170 fGener = generator;
171 Float_t c1min, c1max, c2min, c2max;
172 Int_t n1 = fGener->NCoor1();
173 Int_t n2 = fGener->NCoor2();
174
175 fGener->Coor1Range(c1min, c1max);
176 fGener->Coor2Range(c2min, c2max);
177
178 fHistRadl = new TH2F("hradl","Radiation length map",
179 n2, c2min, c2max, n1, c1min, c1max);
180 fHistAbso = new TH2F("habso","Interaction length map",
181 n2, c2min, c2max, n1, c1min, c1max);
182 fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
183 n2, c2min, c2max, n1, c1min, c1max);
184 //
185 //
b13db077 186
e2afb3b6 187 fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
188 fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
189 fDebug = gAlice->GetDebug();
fe4da5cc 190}
191
e2afb3b6 192//_______________________________________________________________________
fe4da5cc 193AliLego::~AliLego()
194{
8918e700 195 //
196 // Destructor
197 //
198 delete fHistRadl;
199 delete fHistAbso;
200 delete fHistGcm2;
8918e700 201 delete fGener;
e2afb3b6 202 delete fVolumesFwd;
203 delete fVolumesBwd;
fe4da5cc 204}
205
e2afb3b6 206//_______________________________________________________________________
dffd31ef 207void AliLego::BeginEvent()
b13db077 208{
8918e700 209 //
210 // --- Set to 0 radiation length, absorption length and g/cm2 ---
211 //
dffd31ef 212 fTotRadl = 0;
213 fTotAbso = 0;
214 fTotGcm2 = 0;
19237f5f 215
a74d2699 216 if (fDebug) {
e2afb3b6 217 if (fErrorCondition) DumpVolumes();
218 fVolumesFwd->Delete();
219 fVolumesBwd->Delete();
220 fStepsForward = 0;
221 fStepsBackward = 0;
222 fErrorCondition = 0;
223 if (gAlice->CurrentTrack() == 0) fStepBack = 0;
a74d2699 224 }
dffd31ef 225}
226
e2afb3b6 227//_______________________________________________________________________
dffd31ef 228void AliLego::FinishEvent()
229{
8918e700 230 //
231 // Finish the event and update the histos
232 //
0b27b8c7 233 Double_t c1, c2;
234 c1 = fGener->CurCoor1();
235 c2 = fGener->CurCoor2();
10bd7535 236 fHistRadl->Fill(c2,c1,fTotRadl);
237 fHistAbso->Fill(c2,c1,fTotAbso);
238 fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077 239}
240
e2afb3b6 241//_______________________________________________________________________
dffd31ef 242void AliLego::FinishRun()
243{
8918e700 244 //
245 // Store histograms in current Root file
246 //
dffd31ef 247 fHistRadl->Write();
248 fHistAbso->Write();
249 fHistGcm2->Write();
dffd31ef 250
251 // Delete histograms from memory
252 fHistRadl->Delete(); fHistRadl=0;
253 fHistAbso->Delete(); fHistAbso=0;
254 fHistGcm2->Delete(); fHistGcm2=0;
a74d2699 255 //
256 if (fErrorCondition) DumpVolumes();
dffd31ef 257}
258
e2afb3b6 259//_______________________________________________________________________
260void AliLego::Copy(AliLego&) const
8918e700 261{
262 //
263 // Copy *this onto lego -- not implemented
264 //
265 Fatal("Copy","Not implemented!\n");
266}
dffd31ef 267
e2afb3b6 268//_______________________________________________________________________
269void AliLego::StepManager()
270{
271 //
272 // called from AliRun::Stepmanager from gustep.
273 // Accumulate the 3 parameters step by step
274 //
275 static Float_t t;
276 Float_t a,z,dens,radl,absl;
277 Int_t i, id, copy;
278 const char* vol;
279 static Float_t vect[3], dir[3];
280
281 TString tmp1, tmp2;
282 copy = 1;
283 id = gMC->CurrentVolID(copy);
284 vol = gMC->VolName(id);
285 Float_t step = gMC->TrackStep();
286
287 TLorentzVector pos, mom;
288 gMC->TrackPosition(pos);
289 gMC->TrackMomentum(mom);
290
291 Int_t status = 0;
292 if (gMC->IsTrackEntering()) status = 1;
293 if (gMC->IsTrackExiting()) status = 2;
294
7eb618fa 295 if (! fStepBack) {
e2afb3b6 296 // printf("\n volume %s %d", vol, status);
297 //
298 // Normal Forward stepping
299 //
300 if (fDebug) {
301 // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
302
303 //
304 // store volume if different from previous
305 //
7eb618fa 306
a74d2699 307 TClonesArray &lvols = *fVolumesFwd;
7eb618fa 308 if (fStepsForward > 0) {
e2afb3b6 309 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
310 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
a74d2699 311 fStepsForward -= 2;
5f3dc1c9 312 fVolumesFwd->RemoveAt(fStepsForward);
313 fVolumesFwd->RemoveAt(fStepsForward+1);
e2afb3b6 314 }
7eb618fa 315 }
e2afb3b6 316
a74d2699 317 new(lvols[fStepsForward++])
e2afb3b6 318 AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
7eb618fa 319
e2afb3b6 320 } // Debug
321 //
322 // Get current material properties
323
324 gMC->CurrentMaterial(a,z,dens,radl,absl);
325
326 if (z < 1) return;
327
328 // --- See if we have to stop now
329 if (TMath::Abs(pos[2]) > fGener->ZMax() ||
330 pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
331 if (!gMC->IsNewTrack()) {
332 // Not the first step, add past contribution
333 fTotAbso += t/absl;
334 fTotRadl += t/radl;
335 fTotGcm2 += t*dens;
336 if (fDebug) {
337 //
338 // generate "mirror" particle flying back
339 //
340 fStepsBackward = fStepsForward;
341
342 Float_t pmom[3], orig[3];
343 Float_t polar[3] = {0.,0.,0.};
344 Int_t ntr;
345 pmom[0] = -dir[0];
346 pmom[1] = -dir[1];
347 pmom[2] = -dir[2];
348 orig[0] = vect[0];
349 orig[1] = vect[1];
350 orig[2] = vect[2];
351
352 gAlice->SetTrack(1, gAlice->CurrentTrack(),
353 0, pmom, orig, polar, 0., kPNoProcess, ntr);
354 } // debug
355
356 } // not a new track !
7eb618fa 357
e2afb3b6 358 if (fDebug) fStepBack = 1;
359 gMC->StopTrack();
360 return;
361 } // outside scoring region ?
362
363 // --- See how long we have to go
364 for(i=0;i<3;++i) {
365 vect[i]=pos[i];
366 dir[i]=mom[i];
367 }
368
369 t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
370
371 if(step) {
372
373 fTotAbso += step/absl;
374 fTotRadl += step/radl;
375 fTotGcm2 += step*dens;
376 }
377
7eb618fa 378 } else {
379 if (fDebug) {
e2afb3b6 380 //
381 // Geometry debugging
382 // Fly back and compare volume sequence
383 //
384 TClonesArray &lvols = *fVolumesBwd;
385 if (fStepsBackward < fStepsForward) {
386 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
387 if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
388 fStepsBackward += 2;
389 fVolumesBwd->RemoveAt(fStepsBackward-1);
390 fVolumesBwd->RemoveAt(fStepsBackward-2);
7eb618fa 391 }
e2afb3b6 392 }
393
394 fStepsBackward--;
395 // printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
396 if (fStepsBackward < 0) {
7eb618fa 397 gMC->StopTrack();
a74d2699 398 fStepBack = 0;
7eb618fa 399 return;
e2afb3b6 400 }
401
402 new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
403
404 AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
405 if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
406 {
407 printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
408 fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
409 fErrorCondition = 1;
410 }
7eb618fa 411 } // Debug
412 } // bwd/fwd
b13db077 413}
414
e2afb3b6 415//_______________________________________________________________________
a74d2699 416void AliLego::DumpVolumes()
417{
e2afb3b6 418 //
419 // Dump volume sequence in case of error
420 //
421 printf("\n Dumping Volume Sequence:");
422 printf("\n ==============================================");
423
424 for (Int_t i = fStepsForward-1; i>=0; i--)
a74d2699 425 {
e2afb3b6 426 AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
427 AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
428 if (tmp1)
429 printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
430 , i,
431 tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
432 tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
433 if (tmp2 && i>= fStepsBackward)
434 printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
435 , i,
436 tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
437 tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
438
439 printf("\n ............................................................................\n");
a74d2699 440 }
e2afb3b6 441 printf("\n ==============================================\n");
a74d2699 442}