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