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