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