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