AliRoot Core  3dc7879 (3dc7879)
AliLego.cxx
Go to the documentation of this file.
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 /* $Id$ */
17 
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 /*
42 <img src="picts/alilego.gif">
43 */
44 //End_Html
45 //
47 
48 #include "TClonesArray.h"
49 #include "TH2.h"
50 #include "TMath.h"
51 #include "TString.h"
52 #include "TVirtualMC.h"
53 
54 #include "AliLog.h"
55 #include "AliDebugVolume.h"
56 #include "AliLego.h"
57 #include "AliLegoGenerator.h"
58 #include "AliMC.h"
59 #include "AliRun.h"
60 
61 ClassImp(AliLego)
62 
63 //_______________________________________________________________________
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  fRZR(0),
74  fRZA(0),
75  fRZG(0),
76  fVolumesFwd(0),
77  fVolumesBwd(0),
78  fStepBack(0),
79  fStepsBackward(0),
80  fStepsForward(0),
81  fErrorCondition(0),
82  fDebug(0),
83  fStopped(0)
84 {
85  //
86  // Default constructor
87  //
88 }
89 
90 //_______________________________________________________________________
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),
101  fRZR(0),
102  fRZA(0),
103  fRZG(0),
104  fVolumesFwd(0),
105  fVolumesBwd(0),
106  fStepBack(0),
107  fStepsBackward(0),
108  fStepsForward(0),
109  fErrorCondition(0),
110  fDebug(0),
111  fStopped(0)
112 {
113  //
114  // Copy constructor
115  //
116  lego.Copy(*this);
117 }
118 
119 
120 //_______________________________________________________________________
121 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
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),
133  fRZR(0),
134  fRZA(0),
135  fRZG(0),
136  fVolumesFwd(0),
137  fVolumesBwd(0),
138  fStepBack(0),
139  fStepsBackward(0),
140  fStepsForward(0),
141  fErrorCondition(0),
142  fDebug(0),
143  fStopped(0)
144 {
145  //
146  // specify the angular limits and the size of the rectangular box
147  //
148  fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
149  nphi, phimin, phimax, rmin, rmax, zmax);
150  fHistRadl = new TH2F("hradl","Radiation length map",
151  ntheta,thetamin,thetamax,nphi,phimin,phimax);
152  fHistAbso = new TH2F("habso","Interaction length map",
153  ntheta,thetamin,thetamax,nphi,phimin,phimax);
154  fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
155  ntheta,thetamin,thetamax,nphi,phimin,phimax);
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 
165 //
166  fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
167  fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
168  fDebug = AliDebugLevel();
169 }
170 
171 //_______________________________________________________________________
172 AliLego::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),
182  fRZR(0),
183  fRZA(0),
184  fRZG(0),
185  fVolumesFwd(0),
186  fVolumesBwd(0),
187  fStepBack(0),
188  fStepsBackward(0),
189  fStepsForward(0),
190  fErrorCondition(0),
191  fDebug(0),
192  fStopped(0)
193 {
194  //
195  // specify the angular limits and the size of the rectangular box
196  //
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);
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)");
222  //
223  //
224 
225  fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
226  fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
227  fDebug = AliDebugLevel();
228 }
229 
230 //_______________________________________________________________________
232 {
233  //
234  // Destructor
235  //
236  delete fHistRadl;
237  delete fHistAbso;
238  delete fHistGcm2;
239  delete fGener;
240  delete fVolumesFwd;
241  delete fVolumesBwd;
242 }
243 
244 //_______________________________________________________________________
246 {
247  //
248  // --- Set to 0 radiation length, absorption length and g/cm2 ---
249  //
250  fTotRadl = 0;
251  fTotAbso = 0;
252  fTotGcm2 = 0;
253  fStopped = 0;
254 
255  if (fDebug) {
257  fVolumesFwd->Delete();
258  fVolumesBwd->Delete();
259  fStepsForward = 0;
260  fStepsBackward = 0;
261  fErrorCondition = 0;
262  if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
263  }
264 }
265 
266 //_______________________________________________________________________
268 {
269  //
270  // Finish the event and update the histos
271  //
272  Double_t c1, c2;
273  c1 = fGener->CurCoor1();
274  c2 = fGener->CurCoor2();
275  fHistRadl->Fill(c2,c1,fTotRadl);
276  fHistAbso->Fill(c2,c1,fTotAbso);
277  fHistGcm2->Fill(c2,c1,fTotGcm2);
278 }
279 
280 //_______________________________________________________________________
282 {
283  //
284  // Store histograms in current Root file
285  //
286  fHistRadl->Write();
287  fHistAbso->Write();
288  fHistGcm2->Write();
289  fRZR->Write();
290  fRZA->Write();
291  fRZG->Write();
292 
293  // Delete histograms from memory
294  fHistRadl->Delete(); fHistRadl=0;
295  fHistAbso->Delete(); fHistAbso=0;
296  fHistGcm2->Delete(); fHistGcm2=0;
297  //
299 }
300 
301 //_______________________________________________________________________
302 void AliLego::Copy(TObject&) const
303 {
304  //
305  // Copy *this onto lego -- not implemented
306  //
307  AliFatal("Not implemented!");
308 }
309 
310 //_______________________________________________________________________
312 {
313  //
314  // called from AliRun::Stepmanager from gustep.
315  // Accumulate the 3 parameters step by step
316  //
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;
325  id = TVirtualMC::GetMC()->CurrentVolID(copy);
326  vol = TVirtualMC::GetMC()->VolName(id);
327  Float_t step = TVirtualMC::GetMC()->TrackStep();
328 
329  TLorentzVector pos, mom;
330  TVirtualMC::GetMC()->TrackPosition(pos);
331  TVirtualMC::GetMC()->TrackMomentum(mom);
332 
333  Int_t status = 0;
334  if (TVirtualMC::GetMC()->IsTrackEntering()) status = 1;
335  if (TVirtualMC::GetMC()->IsTrackExiting()) status = 2;
336 
337  if (! fStepBack) {
338  // printf("\n volume %s %d", vol, status);
339  //
340  // Normal Forward stepping
341  //
342  if (fDebug) {
343 // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
344 
345  //
346  // store volume if different from previous
347  //
348 
349  TClonesArray &lvols = *fVolumesFwd;
350  if (fStepsForward > 0) {
351  AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
352  if (tmp && tmp->IsVEqual(vol, copy) && TVirtualMC::GetMC()->IsTrackEntering()) {
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 
366  TVirtualMC::GetMC()->CurrentMaterial(a,z,dens,radl,absl);
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 
372 
373  if (z < 1) return;
374 
375  // --- See if we have to stop now
376  if (TMath::Abs(pos[2]) > TMath::Abs(fGener->ZMax()) ||
377  r > fGener->RadMax()) {
378  if (!TVirtualMC::GetMC()->IsNewTrack()) {
379  // Not the first step, add past contribution
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",
388 // pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, TVirtualMC::GetMC()->CurrentVolName(), fTotRadl);
389  if (fDebug) {
390  //
391  // generate "mirror" particle flying back
392  //
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 
406  0, pmom, orig, polar, 0., kPNoProcess, ntr);
407  } // debug
408 
409  } // not a new track !
410 
411  if (fDebug) fStepBack = 1;
412  fStopped = kTRUE;
413  TVirtualMC::GetMC()->StopTrack();
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;
430 // printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n",
431 // pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, TVirtualMC::GetMC()->CurrentVolName(), fTotRadl);
432  }
433 
434  } else {
435  if (fDebug) {
436  //
437  // Geometry debugging
438  // Fly back and compare volume sequence
439  //
440  TClonesArray &lvols = *fVolumesBwd;
442  AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
443  if (tmp && tmp->IsVEqual(vol, copy) && TVirtualMC::GetMC()->IsTrackEntering()) {
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) {
453  TVirtualMC::GetMC()->StopTrack();
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]);
461  if (tmp && !(tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
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
469 }
470 
471 //_______________________________________________________________________
473 {
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--)
481  {
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");
496  }
497  printf("\n ==============================================\n");
498 }
void DumpVolumes()
Definition: AliLego.cxx:472
Int_t fDebug
Error condition flag.
Definition: AliLego.h:60
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TH2F * fHistGcm2
Definition: AliLego.h:49
TClonesArray * fVolumesBwd
Volume sequence forward.
Definition: AliLego.h:55
#define AliDebugLevel()
Float_t fTotRadl
Definition: AliLego.h:44
Bool_t IsVEqual(const char *name, Int_t copy) const
virtual void Coor2Range(Float_t &c2min, Float_t &c2max) const
Float_t fTotAbso
Total Radiation length.
Definition: AliLego.h:45
TH2F * fRZG
Definition: AliLego.h:53
AliLegoGenerator * fGener
Definition: AliLego.h:43
virtual Float_t ZMax() const
TH2F * fHistReta
Definition: AliLego.h:50
virtual void FinishEvent()
Definition: AliLego.cxx:267
Float_t fTotGcm2
Total absorption length.
Definition: AliLego.h:46
virtual void FinishRun()
Definition: AliLego.cxx:281
virtual Int_t NCoor2() const
Int_t fStepBack
Volume sequence backward.
Definition: AliLego.h:56
#define AliWarning(message)
Definition: AliLog.h:541
#define ToAliDebug(level, whatever)
Definition: AliLog.h:379
Float_t PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z)
Float_t Y() const
virtual ~AliLego()
Definition: AliLego.cxx:231
virtual Float_t RadMax() const
Bool_t fStopped
Definition: AliLego.h:61
#define ToAliError(whatever)
Definition: AliLog.h:621
virtual Int_t NCoor1() const
virtual void BeginEvent()
Definition: AliLego.cxx:245
TH2F * fHistRadl
Total G/CM2 traversed.
Definition: AliLego.h:47
Float_t X() const
void polar()
Definition: polar.C:2
Int_t fErrorCondition
Counts steps backward.
Definition: AliLego.h:59
void Copy(TObject &lego) const
Definition: AliLego.cxx:302
AliRun * gAlice
Definition: AliRun.cxx:62
TH2F * fRZR
Definition: AliLego.h:51
Int_t GetCurrentTrackNumber() const
Definition: AliMC.cxx:1846
#define AliFatal(message)
Definition: AliLog.h:640
virtual void PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom, const Float_t *vpos, const Float_t *polar, Float_t tof, TMCProcess mech, Int_t &ntr, Float_t weight=1, Int_t is=0) const
Definition: AliMC.cxx:1925
Int_t CopyNumber() const
TH2F * fHistAbso
Definition: AliLego.h:48
Float_t Step() const
AliLego()
Definition: AliLego.cxx:64
Int_t fStepsForward
Counts steps forward.
Definition: AliLego.h:58
TClonesArray * fVolumesFwd
Definition: AliLego.h:54
virtual void Coor1Range(Float_t &c1min, Float_t &c1max) const
const char * Status() const
TH2F * fRZA
Definition: AliLego.h:52
virtual void StepManager()
Definition: AliLego.cxx:311
Int_t fStepsBackward
Flag for backstepping.
Definition: AliLego.h:57
virtual Float_t CurCoor2() const
AliMC * GetMCApp() const
Definition: AliRun.h:53
Float_t Z() const
virtual Float_t CurCoor1() const