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