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