prevent running if CDB snapshot setting failed
[u/mrichter/AliRoot.git] / STEER / 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   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 //_______________________________________________________________________
91 AliLego::AliLego(const AliLego &lego):
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 //_______________________________________________________________________
231 AliLego::~AliLego()
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 //_______________________________________________________________________
245 void AliLego::BeginEvent()
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) {
256     if (fErrorCondition) ToAliDebug(1, DumpVolumes());
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 //_______________________________________________________________________
267 void AliLego::FinishEvent()
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 //_______________________________________________________________________
281 void AliLego::FinishRun()
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   //
298   if (fErrorCondition) ToAliError(DumpVolumes());
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 //_______________________________________________________________________
311 void AliLego::StepManager() 
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                     //
393                     fStepsBackward = fStepsForward;
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                     
405                     gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(), 
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;
441             if (fStepsBackward < fStepsForward) {
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 //_______________________________________________________________________
472 void AliLego::DumpVolumes()
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 }