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