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