Use fStepBack = 1 only in debug mode.
[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 /*
17 $Log$
18 Revision 1.23  2001/05/20 10:10:39  morsch
19 - Debug output at the beginning of new event and end of run.
20 - Filter out boundary loops.
21
22 Revision 1.22  2001/05/11 13:22:40  morsch
23 If run with debug option (from gAlice) geantinos are sent back and volume sequence forward/backward is compared.
24 Can be very verbous in some cases.
25
26 Revision 1.21  2000/12/15 10:33:59  morsch
27 Invert coordinates to make meaningful zylindrical plots.
28
29 Revision 1.20  2000/11/30 07:12:48  alibrary
30 Introducing new Rndm and QA classes
31
32 Revision 1.19  2000/10/26 14:13:05  morsch
33 - Change from coordinates theta, phi to general coordinates Coor1 and Coor2.
34 - Lego generator instance can be passed in constructor.
35
36 Revision 1.18  2000/10/02 21:28:14  fca
37 Removal of useless dependecies via forward declarations
38
39 Revision 1.17  2000/07/12 08:56:25  fca
40 Coding convention correction and warning removal
41
42 Revision 1.16  2000/05/26 08:35:03  fca
43 Move the check on z after z has been retrieved
44
45 Revision 1.15  2000/05/16 13:10:40  fca
46 New method IsNewTrack and fix for a problem in Father-Daughter relations
47
48 Revision 1.14  2000/04/27 10:38:21  fca
49 Correct termination of Lego Run and introduce Lego getter in AliRun
50
51 Revision 1.13  2000/04/26 10:17:31  fca
52 Changes in Lego for G4 compatibility
53
54 Revision 1.12  2000/04/07 11:12:33  fca
55 G4 compatibility changes
56
57 Revision 1.11  2000/03/22 13:42:26  fca
58 SetGenerator does not replace an existing generator, ResetGenerator does
59
60 Revision 1.10  2000/02/23 16:25:22  fca
61 AliVMC and AliGeant3 classes introduced
62 ReadEuclid moved from AliRun to AliModule
63
64 Revision 1.9  1999/12/03 10:54:01  fca
65 Fix lego summary
66
67 Revision 1.8  1999/10/01 09:54:33  fca
68 Correct logics for Lego StepManager
69
70 Revision 1.7  1999/09/29 09:24:29  fca
71 Introduction of the Copyright and cvs Log
72
73 */
74
75 //////////////////////////////////////////////////////////////
76 //////////////////////////////////////////////////////////////
77 //                            
78 //  Utility class to evaluate the material budget from
79 //  a given radius to the surface of an arbitrary cylinder
80 //  along radial directions from the centre:
81 // 
82 //   - radiation length
83 //   - Interaction length
84 //   - g/cm2
85 // 
86 //  Geantinos are shot in the bins in the fNtheta bins in theta
87 //  and fNphi bins in phi with specified rectangular limits.
88 //  The statistics are accumulated per
89 //    fRadMin < r < fRadMax    and  <0 < z < fZMax
90 //
91 //  To activate this option, you can do:
92 //      Root > gAlice.RunLego();
93 //  or  Root > .x menu.C  then select menu item "RunLego"
94 //  Note that when calling gAlice->RunLego, an optional list
95 //  of arguments may be specified.
96 //
97 //Begin_Html
98 /*
99 <img src="picts/alilego.gif">
100 */
101 //End_Html
102 //
103 //////////////////////////////////////////////////////////////
104
105 #include "TMath.h"
106
107 #include "AliLego.h"
108
109 #include "AliDebugVolume.h"
110 #include "AliRun.h"
111 #include "AliLegoGenerator.h"
112 #include "AliConst.h"
113 #include "AliMC.h"
114 #include "TH2.h"
115 #include "../TGeant3/TGeant3.h"
116 #include "TString.h"
117 #include "TClonesArray.h"
118
119 ClassImp(AliLego)
120
121
122 //___________________________________________
123 AliLego::AliLego()
124 {
125   //
126   // Default constructor
127   //
128   fHistRadl = 0;
129   fHistAbso = 0;
130   fHistGcm2 = 0;
131   fHistReta = 0;
132 }
133
134 //___________________________________________
135 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, 
136                  Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
137                  Float_t rmin, Float_t rmax, Float_t zmax)
138   : TNamed("Lego Generator",title)
139 {
140   //
141   // specify the angular limits and the size of the rectangular box
142   //
143    fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
144                        nphi, phimin, phimax, rmin, rmax, zmax);
145    
146
147    
148    fHistRadl = new TH2F("hradl","Radiation length map",    
149                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
150    fHistAbso = new TH2F("habso","Interaction length map",  
151                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
152    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
153                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
154 //
155    fStepBack       = 0;
156    fStepsBackward  = 0;
157    fStepsForward   = 0;
158    fStepBack       = 0;
159    fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
160    fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
161    fDebug          = gAlice->GetDebug();
162    fErrorCondition = 0;   
163 }
164
165 AliLego::AliLego(const char *title, AliLegoGenerator* generator)
166   : TNamed("Lego Generator",title)
167 {
168   //
169   // specify the angular limits and the size of the rectangular box
170   //
171    fGener = generator;
172    Float_t c1min, c1max, c2min, c2max;
173    Int_t n1 = fGener->NCoor1();
174    Int_t n2 = fGener->NCoor2();
175    
176    fGener->Coor1Range(c1min, c1max);
177    fGener->Coor2Range(c2min, c2max);   
178
179    fHistRadl = new TH2F("hradl","Radiation length map",    
180                         n2, c2min, c2max, n1, c1min, c1max);
181    fHistAbso = new TH2F("habso","Interaction length map",  
182                         n2, c2min, c2max, n1, c1min, c1max);
183    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
184                         n2, c2min, c2max, n1, c1min, c1max);
185 //
186 //
187    fStepBack       = 0;
188    fStepsBackward  = 0;
189    fStepsForward   = 0;
190    fStepBack       = 0;
191    fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
192    fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
193    fDebug          = gAlice->GetDebug();
194    fErrorCondition =0;
195 }
196
197 //___________________________________________
198 AliLego::~AliLego()
199 {
200   //
201   // Destructor
202   //
203   delete fHistRadl;
204   delete fHistAbso;
205   delete fHistGcm2;
206   delete fGener;
207   if (fVolumesFwd) delete fVolumesFwd;
208   if (fVolumesBwd) delete fVolumesBwd;
209 }
210
211 //___________________________________________
212 void AliLego::BeginEvent()
213 {
214   //
215   // --- Set to 0 radiation length, absorption length and g/cm2 ---
216   //
217   fTotRadl = 0;
218   fTotAbso = 0;
219   fTotGcm2 = 0;
220
221   if (fDebug) {
222       if (fErrorCondition) DumpVolumes();
223       fVolumesFwd->Delete();
224       fVolumesBwd->Delete();
225       fStepsForward    = 0;
226       fStepsBackward   = 0;               
227       fErrorCondition  = 0;
228       if (gAlice->CurrentTrack() == 0) fStepBack = 0;
229   }
230 }
231
232 //___________________________________________
233 void AliLego::FinishEvent()
234 {
235   //
236   // Finish the event and update the histos
237   //
238   Double_t c1, c2;
239   c1 = fGener->CurCoor1();
240   c2 = fGener->CurCoor2();
241   fHistRadl->Fill(c2,c1,fTotRadl);
242   fHistAbso->Fill(c2,c1,fTotAbso);
243   fHistGcm2->Fill(c2,c1,fTotGcm2);
244 }
245
246 //___________________________________________
247 void AliLego::FinishRun()
248 {
249   //
250   // Store histograms in current Root file
251   //
252   fHistRadl->Write();
253   fHistAbso->Write();
254   fHistGcm2->Write();
255
256   // Delete histograms from memory
257   fHistRadl->Delete(); fHistRadl=0;
258   fHistAbso->Delete(); fHistAbso=0;
259   fHistGcm2->Delete(); fHistGcm2=0;
260   //
261   if (fErrorCondition) DumpVolumes();
262 }
263
264 //___________________________________________
265 void AliLego::Copy(AliLego &lego) const
266 {
267   //
268   // Copy *this onto lego -- not implemented
269   //
270   Fatal("Copy","Not implemented!\n");
271 }
272
273 //___________________________________________
274 void AliLego::StepManager() {
275 // called from AliRun::Stepmanager from gustep.
276 // Accumulate the 3 parameters step by step
277
278     static Float_t t;
279     Float_t a,z,dens,radl,absl;
280     Int_t i, id, copy;
281     const char* vol;
282     static Float_t vect[3], dir[3];
283     
284     TString tmp1, tmp2;
285     copy = 1;
286     id  = gMC->CurrentVolID(copy);
287     vol = gMC->VolName(id);
288     Float_t step  = gMC->TrackStep();
289
290    TLorentzVector pos, mom; 
291    gMC->TrackPosition(pos);  
292    gMC->TrackMomentum(mom);
293    
294    Int_t status = 0;
295    if (gMC->IsTrackEntering()) status = 1;
296    if (gMC->IsTrackExiting())  status = 2;
297
298    
299
300    
301   if (! fStepBack) {
302 //      printf("\n volume %s %d", vol, status);      
303 //
304 // Normal Forward stepping
305 //
306       if (fDebug) {
307 //        printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);          
308
309 //
310 // store volume if different from previous
311 //
312           
313           TClonesArray &lvols = *fVolumesFwd;
314           if (fStepsForward > 0) {
315               AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesFwd)[fStepsForward-1];         
316               if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) {
317                   step += ((AliDebugVolume*) lvols[fStepsForward])->Step();
318                   fStepsForward -= 2;
319                   delete ((AliDebugVolume*) lvols[fStepsForward]);
320                   delete ((AliDebugVolume*) lvols[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);
332       
333       if (z < 1) return;
334       
335 // --- See if we have to stop now
336        if (TMath::Abs(pos[2]) > fGener->ZMax()  || 
337            pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
338            if (!gMC->IsNewTrack()) {
339                // Not the first step, add past contribution
340                fTotAbso += t/absl;
341                fTotRadl += t/radl;
342                fTotGcm2 += t*dens;
343                if (fDebug) {
344 //
345 //  generate "mirror" particle flying back
346 //
347                    fStepsBackward = fStepsForward;
348                
349                    Float_t pmom[3], orig[3];
350                    Float_t polar[3] = {0.,0.,0.};
351                    Int_t ntr;
352                    pmom[0] = -dir[0];
353                    pmom[1] = -dir[1];      
354                    pmom[2] = -dir[2];
355                    orig[0] =  vect[0];
356                    orig[1] =  vect[1];     
357                    orig[2] =  vect[2];
358                    
359                    gAlice->SetTrack(1, gAlice->CurrentTrack(), 
360                                     0, pmom, orig, polar, 0., kPNoProcess, ntr);
361                } // debug
362                
363            } // not a new track !
364
365            if (fDebug) fStepBack = 1;
366            gMC->StopTrack();
367            return;
368        } // outside scoring region ?
369        
370 // --- See how long we have to go
371        for(i=0;i<3;++i) {
372            vect[i]=pos[i];
373            dir[i]=mom[i];
374        }
375        
376        t  = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
377        
378        if(step) {
379
380            fTotAbso += step/absl;
381            fTotRadl += step/radl;
382            fTotGcm2 += step*dens;
383        }
384        
385
386   } else {
387       if (fDebug) {
388 //
389 // Geometry debugging
390 // Fly back and compare volume sequence
391 //
392           TClonesArray &lvols = *fVolumesBwd;
393           if (fStepsBackward < fStepsForward) {
394               AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesBwd)[fStepsBackward];          
395               if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) {
396                   fStepsBackward += 2;
397                   delete ((AliDebugVolume*) lvols[fStepsBackward-1]);
398                   delete ((AliDebugVolume*) lvols[fStepsBackward-2]);
399                   step += ((AliDebugVolume*) lvols[fStepsBackward])->Step();
400               }
401           } 
402
403           fStepsBackward--;
404 //        printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);      
405           if (fStepsBackward < 0) {
406               gMC->StopTrack();
407               fStepBack = 0;
408               return;
409           }
410           
411           new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
412           
413           AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesFwd)[fStepsBackward];
414           if (! (tmp->IsEqual(vol, copy)) && (!fErrorCondition)) 
415           {
416               printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n", 
417                      fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
418               fErrorCondition = 1;
419           } 
420       } // Debug
421   } // bwd/fwd
422 }
423
424 void AliLego::DumpVolumes()
425 {
426 //
427 // Dump volume sequence in case of error
428 //
429     printf("\n Dumping Volume Sequence:");
430     printf("\n ==============================================");
431     
432     for (Int_t i = fStepsForward-1; i>=0; i--)
433     {
434         AliDebugVolume* tmp1 = (AliDebugVolume*) (*fVolumesFwd)[i];
435         AliDebugVolume* tmp2 = (AliDebugVolume*) (*fVolumesBwd)[i];
436         printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s\n"
437                , i, 
438                tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(), 
439                tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
440         
441         printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s\n"
442                , i, 
443                tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(), 
444                tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
445         
446         printf("\n ............................................................................\n");
447     }
448     printf("\n ==============================================\n");
449 }
450
451
452
453
454
455
456