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