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