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 /* $Id$ */
  17
  18 //////////////////////////////////////////////////////////////
  19 //////////////////////////////////////////////////////////////
  20 //
  21 //  Utility class to evaluate the material budget from
  22 //  a given radius to the surface of an arbitrary cylinder
  23 //  along radial directions from the centre:
  24 //
  25 //   - radiation length
  26 //   - Interaction length
  27 //   - g/cm2
  28 //
  29 //  Geantinos are shot in the bins in the fNtheta bins in theta
  30 //  and fNphi bins in phi with specified rectangular limits.
  31 //  The statistics are accumulated per
  32 //    fRadMin < r < fRadMax    and  <0 < z < fZMax
  33 //
  34 //  To activate this option, you can do:
  35 //      Root > gAlice.RunLego();
  36 //  or  Root > .x menu.C  then select menu item "RunLego"
  37 //  Note that when calling gAlice->RunLego, an optional list
  38 //  of arguments may be specified.
  39 //
  40 //Begin_Html
  41 /*
  42 <img src="picts/alilego.gif">
  43 */
  44 //End_Html
  45 //
  46 //////////////////////////////////////////////////////////////
  47
  48 #include "TClonesArray.h"
  49 #include "TH2.h"
  50 #include "TMath.h"
  51 #include "TString.h"
  52
  53 #include "AliConst.h"
  54 #include "AliDebugVolume.h"
  55 #include "AliLego.h"
  56 #include "AliLegoGenerator.h"
  57 #include "AliRun.h"
  58
  59 ClassImp(AliLego)
  60
  61 //_______________________________________________________________________
  62 AliLego::AliLego():
  63   fGener(0),
  64   fTotRadl(0),
  65   fTotAbso(0),
  66   fTotGcm2(0),
  67   fHistRadl(0),
  68   fHistAbso(0),
  69   fHistGcm2(0),
  70   fHistReta(0),
  71   fVolumesFwd(0),
  72   fVolumesBwd(0),
  73   fStepBack(0),
  74   fStepsBackward(0),
  75   fStepsForward(0),
  76   fErrorCondition(0),
  77   fDebug(0)
  78 {
  79   //
  80   // Default constructor
  81   //
  82 }
  83
  84 //_______________________________________________________________________
  85 AliLego::AliLego(const AliLego &lego):
  86   TNamed(lego),
  87   fGener(0),
  88   fTotRadl(0),
  89   fTotAbso(0),
  90   fTotGcm2(0),
  91   fHistRadl(0),
  92   fHistAbso(0),
  93   fHistGcm2(0),
  94   fHistReta(0),
  95   fVolumesFwd(0),
  96   fVolumesBwd(0),
  97   fStepBack(0),
  98   fStepsBackward(0),
  99   fStepsForward(0),
 100   fErrorCondition(0),
 101   fDebug(0)
 102 {
 103   //
 104   // Copy constructor
 105   //
 106   lego.Copy(*this);
 107 }
 108
 109
 110 //_______________________________________________________________________
 111 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
 112                  Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
 113                  Float_t rmin, Float_t rmax, Float_t zmax):
 114   TNamed("Lego Generator",title),
 115   fGener(0),
 116   fTotRadl(0),
 117   fTotAbso(0),
 118   fTotGcm2(0),
 119   fHistRadl(0),
 120   fHistAbso(0),
 121   fHistGcm2(0),
 122   fHistReta(0),
 123   fVolumesFwd(0),
 124   fVolumesBwd(0),
 125   fStepBack(0),
 126   fStepsBackward(0),
 127   fStepsForward(0),
 128   fErrorCondition(0),
 129   fDebug(0)
 130 {
 131   //
 132   // specify the angular limits and the size of the rectangular box
 133   //
 134    fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
 135                        nphi, phimin, phimax, rmin, rmax, zmax);
 136    fHistRadl = new TH2F("hradl","Radiation length map",
 137                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 138    fHistAbso = new TH2F("habso","Interaction length map",
 139                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 140    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
 141                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 142 //
 143    fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
 144    fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
 145    fDebug          = gAlice->GetDebug();
 146 }
 147
 148 //_______________________________________________________________________
 149 AliLego::AliLego(const char *title, AliLegoGenerator* generator):
 150   TNamed("Lego Generator",title),
 151   fGener(0),
 152   fTotRadl(0),
 153   fTotAbso(0),
 154   fTotGcm2(0),
 155   fHistRadl(0),
 156   fHistAbso(0),
 157   fHistGcm2(0),
 158   fHistReta(0),
 159   fVolumesFwd(0),
 160   fVolumesBwd(0),
 161   fStepBack(0),
 162   fStepsBackward(0),
 163   fStepsForward(0),
 164   fErrorCondition(0),
 165   fDebug(0)
 166 {
 167   //
 168   // specify the angular limits and the size of the rectangular box
 169   //
 170   fGener = generator;
 171   Float_t c1min, c1max, c2min, c2max;
 172   Int_t n1 = fGener->NCoor1();
 173   Int_t n2 = fGener->NCoor2();
 174
 175   fGener->Coor1Range(c1min, c1max);
 176   fGener->Coor2Range(c2min, c2max);
 177
 178   fHistRadl = new TH2F("hradl","Radiation length map",
 179                        n2, c2min, c2max, n1, c1min, c1max);
 180   fHistAbso = new TH2F("habso","Interaction length map",
 181                        n2, c2min, c2max, n1, c1min, c1max);
 182   fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
 183                        n2, c2min, c2max, n1, c1min, c1max);
 184   //
 185   //
 186
 187   fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
 188   fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
 189   fDebug          = gAlice->GetDebug();
 190 }
 191
 192 //_______________________________________________________________________
 193 AliLego::~AliLego()
 194 {
 195   //
 196   // Destructor
 197   //
 198   delete fHistRadl;
 199   delete fHistAbso;
 200   delete fHistGcm2;
 201   delete fGener;
 202   delete fVolumesFwd;
 203   delete fVolumesBwd;
 204 }
 205
 206 //_______________________________________________________________________
 207 void AliLego::BeginEvent()
 208 {
 209   //
 210   // --- Set to 0 radiation length, absorption length and g/cm2 ---
 211   //
 212   fTotRadl = 0;
 213   fTotAbso = 0;
 214   fTotGcm2 = 0;
 215
 216   if (fDebug) {
 217     if (fErrorCondition) DumpVolumes();
 218     fVolumesFwd->Delete();
 219     fVolumesBwd->Delete();
 220     fStepsForward    = 0;
 221     fStepsBackward   = 0;
 222     fErrorCondition  = 0;
 223     if (gAlice->GetCurrentTrackNumber() == 0) fStepBack = 0;
 224   }
 225 }
 226
 227 //_______________________________________________________________________
 228 void AliLego::FinishEvent()
 229 {
 230   //
 231   // Finish the event and update the histos
 232   //
 233   Double_t c1, c2;
 234   c1 = fGener->CurCoor1();
 235   c2 = fGener->CurCoor2();
 236   fHistRadl->Fill(c2,c1,fTotRadl);
 237   fHistAbso->Fill(c2,c1,fTotAbso);
 238   fHistGcm2->Fill(c2,c1,fTotGcm2);
 239 }
 240
 241 //_______________________________________________________________________
 242 void AliLego::FinishRun()
 243 {
 244   //
 245   // Store histograms in current Root file
 246   //
 247   fHistRadl->Write();
 248   fHistAbso->Write();
 249   fHistGcm2->Write();
 250
 251   // Delete histograms from memory
 252   fHistRadl->Delete(); fHistRadl=0;
 253   fHistAbso->Delete(); fHistAbso=0;
 254   fHistGcm2->Delete(); fHistGcm2=0;
 255   //
 256   if (fErrorCondition) DumpVolumes();
 257 }
 258
 259 //_______________________________________________________________________
 260 void AliLego::Copy(AliLego&) const
 261 {
 262   //
 263   // Copy *this onto lego -- not implemented
 264   //
 265   Fatal("Copy","Not implemented!\n");
 266 }
 267
 268 //_______________________________________________________________________
 269 void AliLego::StepManager()
 270 {
 271   //
 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   if (! fStepBack) {
 296     //      printf("\n volume %s %d", vol, status);
 297     //
 298     // Normal Forward stepping
 299     //
 300     if (fDebug) {
 301       //          printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
 302
 303       //
 304       // store volume if different from previous
 305       //
 306
 307           TClonesArray &lvols = *fVolumesFwd;
 308           if (fStepsForward > 0) {
 309         AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
 310         if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
 311                   fStepsForward -= 2;
 312                   fVolumesFwd->RemoveAt(fStepsForward);
 313                   fVolumesFwd->RemoveAt(fStepsForward+1);
 314         }
 315           }
 316
 317           new(lvols[fStepsForward++])
 318         AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
 319
 320     } // Debug
 321     //
 322     // Get current material properties
 323
 324     gMC->CurrentMaterial(a,z,dens,radl,absl);
 325
 326     if (z < 1) return;
 327
 328     // --- See if we have to stop now
 329     if (TMath::Abs(pos[2]) > fGener->ZMax()  ||
 330         pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
 331       if (!gMC->IsNewTrack()) {
 332         // Not the first step, add past contribution
 333         fTotAbso += t/absl;
 334         fTotRadl += t/radl;
 335         fTotGcm2 += t*dens;
 336         if (fDebug) {
 337           //
 338           //  generate "mirror" particle flying back
 339           //
 340           fStepsBackward = fStepsForward;
 341
 342           Float_t pmom[3], orig[3];
 343           Float_t polar[3] = {0.,0.,0.};
 344           Int_t ntr;
 345           pmom[0] = -dir[0];
 346           pmom[1] = -dir[1];
 347           pmom[2] = -dir[2];
 348           orig[0] =  vect[0];
 349           orig[1] =  vect[1];
 350           orig[2] =  vect[2];
 351
 352           gAlice->PushTrack(1, gAlice->GetCurrentTrackNumber(),
 353                            0, pmom, orig, polar, 0., kPNoProcess, ntr);
 354         } // debug
 355
 356       } // not a new track !
 357
 358       if (fDebug) fStepBack = 1;
 359       gMC->StopTrack();
 360       return;
 361     } // outside scoring region ?
 362
 363     // --- See how long we have to go
 364     for(i=0;i<3;++i) {
 365       vect[i]=pos[i];
 366       dir[i]=mom[i];
 367     }
 368
 369     t  = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
 370
 371     if(step) {
 372
 373       fTotAbso += step/absl;
 374       fTotRadl += step/radl;
 375       fTotGcm2 += step*dens;
 376     }
 377
 378   } else {
 379       if (fDebug) {
 380         //
 381         // Geometry debugging
 382         // Fly back and compare volume sequence
 383         //
 384         TClonesArray &lvols = *fVolumesBwd;
 385         if (fStepsBackward < fStepsForward) {
 386               AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
 387               if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
 388             fStepsBackward += 2;
 389             fVolumesBwd->RemoveAt(fStepsBackward-1);
 390             fVolumesBwd->RemoveAt(fStepsBackward-2);
 391               }
 392         }
 393
 394         fStepsBackward--;
 395         //        printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
 396         if (fStepsBackward < 0) {
 397               gMC->StopTrack();
 398               fStepBack = 0;
 399               return;
 400         }
 401
 402         new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
 403
 404         AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
 405         if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
 406           {
 407             printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
 408                    fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
 409             fErrorCondition = 1;
 410           }
 411       } // Debug
 412   } // bwd/fwd
 413 }
 414
 415 //_______________________________________________________________________
 416 void AliLego::DumpVolumes()
 417 {
 418   //
 419   // Dump volume sequence in case of error
 420   //
 421   printf("\n Dumping Volume Sequence:");
 422   printf("\n ==============================================");
 423
 424   for (Int_t i = fStepsForward-1; i>=0; i--)
 425     {
 426       AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
 427       AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
 428       if (tmp1)
 429         printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
 430                , i,
 431                tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
 432                tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
 433       if (tmp2 && i>= fStepsBackward)
 434         printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
 435                , i,
 436                tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
 437                tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
 438
 439       printf("\n ............................................................................\n");
 440     }
 441   printf("\n ==============================================\n");
 442 }