Moving to the new VMC naming convention
[u/mrichter/AliRoot.git] / STEER / AliLego.cxx
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
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 /*
42 <img src="picts/alilego.gif">
43 */
44 //End_Html
45 //
46 //////////////////////////////////////////////////////////////
47
48 #include "TClonesArray.h"
49 #include "TH2.h"
50 #include "TMath.h"
51 #include "TString.h"
52
53 #include "AliConst.h"
54 #include "AliDebugVolume.h"
55 #include "AliLego.h"
56 #include "AliLegoGenerator.h"
57 #include "AliRun.h"
58
59 ClassImp(AliLego)
60
61 //_______________________________________________________________________
62 AliLego::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)
78 {
79   //
80   // Default constructor
81   //
82 }
83
84 //_______________________________________________________________________
85 AliLego::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 //_______________________________________________________________________
111 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, 
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)
130 {
131   //
132   // specify the angular limits and the size of the rectangular box
133   //
134    fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
135                        nphi, phimin, phimax, rmin, rmax, zmax);
136    fHistRadl = new TH2F("hradl","Radiation length map",    
137                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
138    fHistAbso = new TH2F("habso","Interaction length map",  
139                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
140    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
141                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
142 //
143    fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
144    fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
145    fDebug          = gAlice->GetDebug();
146 }
147
148 //_______________________________________________________________________
149 AliLego::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)
166 {
167   //
168   // specify the angular limits and the size of the rectangular box
169   //
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   //
186
187   fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
188   fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
189   fDebug          = gAlice->GetDebug();
190 }
191
192 //_______________________________________________________________________
193 AliLego::~AliLego()
194 {
195   //
196   // Destructor
197   //
198   delete fHistRadl;
199   delete fHistAbso;
200   delete fHistGcm2;
201   delete fGener;
202   delete fVolumesFwd;
203   delete fVolumesBwd;
204 }
205
206 //_______________________________________________________________________
207 void AliLego::BeginEvent()
208 {
209   //
210   // --- Set to 0 radiation length, absorption length and g/cm2 ---
211   //
212   fTotRadl = 0;
213   fTotAbso = 0;
214   fTotGcm2 = 0;
215
216   if (fDebug) {
217     if (fErrorCondition) DumpVolumes();
218     fVolumesFwd->Delete();
219     fVolumesBwd->Delete();
220     fStepsForward    = 0;
221     fStepsBackward   = 0;                 
222     fErrorCondition  = 0;
223     if (gAlice->GetCurrentTrackNumber() == 0) fStepBack = 0;
224   }
225 }
226
227 //_______________________________________________________________________
228 void AliLego::FinishEvent()
229 {
230   //
231   // Finish the event and update the histos
232   //
233   Double_t c1, c2;
234   c1 = fGener->CurCoor1();
235   c2 = fGener->CurCoor2();
236   fHistRadl->Fill(c2,c1,fTotRadl);
237   fHistAbso->Fill(c2,c1,fTotAbso);
238   fHistGcm2->Fill(c2,c1,fTotGcm2);
239 }
240
241 //_______________________________________________________________________
242 void AliLego::FinishRun()
243 {
244   //
245   // Store histograms in current Root file
246   //
247   fHistRadl->Write();
248   fHistAbso->Write();
249   fHistGcm2->Write();
250
251   // Delete histograms from memory
252   fHistRadl->Delete(); fHistRadl=0;
253   fHistAbso->Delete(); fHistAbso=0;
254   fHistGcm2->Delete(); fHistGcm2=0;
255   //
256   if (fErrorCondition) DumpVolumes();
257 }
258
259 //_______________________________________________________________________
260 void AliLego::Copy(AliLego&) const
261 {
262   //
263   // Copy *this onto lego -- not implemented
264   //
265   Fatal("Copy","Not implemented!\n");
266 }
267
268 //_______________________________________________________________________
269 void 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   
295   if (! fStepBack) {
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       //
306           
307           TClonesArray &lvols = *fVolumesFwd;
308           if (fStepsForward > 0) {
309         AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
310         if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
311                   fStepsForward -= 2;
312                   fVolumesFwd->RemoveAt(fStepsForward);
313                   fVolumesFwd->RemoveAt(fStepsForward+1);                 
314         }
315           }
316       
317           new(lvols[fStepsForward++]) 
318         AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
319       
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->PushTrack(1, gAlice->GetCurrentTrackNumber(), 
353                            0, pmom, orig, polar, 0., kPNoProcess, ntr);
354         } // debug
355         
356       } // not a new track !
357       
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     
378   } else {
379       if (fDebug) {
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);              
391               }
392         } 
393         
394         fStepsBackward--;
395         //        printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);      
396         if (fStepsBackward < 0) {
397               gMC->StopTrack();
398               fStepBack = 0;
399               return;
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           } 
411       } // Debug
412   } // bwd/fwd
413 }
414
415 //_______________________________________________________________________
416 void AliLego::DumpVolumes()
417 {
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--)
425     {
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");
440     }
441   printf("\n ==============================================\n");
442 }