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