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