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 #include "TVirtualMC.h"
  53
  54 #include "AliDebugVolume.h"
  55 #include "AliLego.h"
  56 #include "AliLegoGenerator.h"
  57 #include "AliMC.h"
  58 #include "AliRun.h"
  59
  60 ClassImp(AliLego)
  61
  62 //_______________________________________________________________________
  63 AliLego::AliLego():
  64   fGener(0),
  65   fTotRadl(0),
  66   fTotAbso(0),
  67   fTotGcm2(0),
  68   fHistRadl(0),
  69   fHistAbso(0),
  70   fHistGcm2(0),
  71   fHistReta(0),
  72   fVolumesFwd(0),
  73   fVolumesBwd(0),
  74   fStepBack(0),
  75   fStepsBackward(0),
  76   fStepsForward(0),
  77   fErrorCondition(0),
  78   fDebug(0)
  79 {
  80   //
  81   // Default constructor
  82   //
  83 }
  84
  85 //_______________________________________________________________________
  86 AliLego::AliLego(const AliLego &lego):
  87   TNamed(lego),
  88   fGener(0),
  89   fTotRadl(0),
  90   fTotAbso(0),
  91   fTotGcm2(0),
  92   fHistRadl(0),
  93   fHistAbso(0),
  94   fHistGcm2(0),
  95   fHistReta(0),
  96   fVolumesFwd(0),
  97   fVolumesBwd(0),
  98   fStepBack(0),
  99   fStepsBackward(0),
 100   fStepsForward(0),
 101   fErrorCondition(0),
 102   fDebug(0)
 103 {
 104   //
 105   // Copy constructor
 106   //
 107   lego.Copy(*this);
 108 }
 109
 110
 111 //_______________________________________________________________________
 112 AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
 113                  Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
 114                  Float_t rmin, Float_t rmax, Float_t zmax):
 115   TNamed("Lego Generator",title),
 116   fGener(0),
 117   fTotRadl(0),
 118   fTotAbso(0),
 119   fTotGcm2(0),
 120   fHistRadl(0),
 121   fHistAbso(0),
 122   fHistGcm2(0),
 123   fHistReta(0),
 124   fVolumesFwd(0),
 125   fVolumesBwd(0),
 126   fStepBack(0),
 127   fStepsBackward(0),
 128   fStepsForward(0),
 129   fErrorCondition(0),
 130   fDebug(0)
 131 {
 132   //
 133   // specify the angular limits and the size of the rectangular box
 134   //
 135    fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
 136                        nphi, phimin, phimax, rmin, rmax, zmax);
 137    fHistRadl = new TH2F("hradl","Radiation length map",
 138                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 139    fHistAbso = new TH2F("habso","Interaction length map",
 140                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 141    fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
 142                         ntheta,thetamin,thetamax,nphi,phimin,phimax);
 143 //
 144    fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
 145    fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
 146    fDebug          = gAlice->GetDebug();
 147 }
 148
 149 //_______________________________________________________________________
 150 AliLego::AliLego(const char *title, AliLegoGenerator* generator):
 151   TNamed("Lego Generator",title),
 152   fGener(0),
 153   fTotRadl(0),
 154   fTotAbso(0),
 155   fTotGcm2(0),
 156   fHistRadl(0),
 157   fHistAbso(0),
 158   fHistGcm2(0),
 159   fHistReta(0),
 160   fVolumesFwd(0),
 161   fVolumesBwd(0),
 162   fStepBack(0),
 163   fStepsBackward(0),
 164   fStepsForward(0),
 165   fErrorCondition(0),
 166   fDebug(0)
 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
 188   fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
 189   fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
 190   fDebug          = gAlice->GetDebug();
 191 }
 192
 193 //_______________________________________________________________________
 194 AliLego::~AliLego()
 195 {
 196   //
 197   // Destructor
 198   //
 199   delete fHistRadl;
 200   delete fHistAbso;
 201   delete fHistGcm2;
 202   delete fGener;
 203   delete fVolumesFwd;
 204   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
 217   if (fDebug) {
 218     if (fErrorCondition) DumpVolumes();
 219     fVolumesFwd->Delete();
 220     fVolumesBwd->Delete();
 221     fStepsForward    = 0;
 222     fStepsBackward   = 0;
 223     fErrorCondition  = 0;
 224     if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
 225   }
 226 }
 227
 228 //_______________________________________________________________________
 229 void AliLego::FinishEvent()
 230 {
 231   //
 232   // Finish the event and update the histos
 233   //
 234   Double_t c1, c2;
 235   c1 = fGener->CurCoor1();
 236   c2 = fGener->CurCoor2();
 237   fHistRadl->Fill(c2,c1,fTotRadl);
 238   fHistAbso->Fill(c2,c1,fTotAbso);
 239   fHistGcm2->Fill(c2,c1,fTotGcm2);
 240 }
 241
 242 //_______________________________________________________________________
 243 void AliLego::FinishRun()
 244 {
 245   //
 246   // Store histograms in current Root file
 247   //
 248   fHistRadl->Write();
 249   fHistAbso->Write();
 250   fHistGcm2->Write();
 251
 252   // Delete histograms from memory
 253   fHistRadl->Delete(); fHistRadl=0;
 254   fHistAbso->Delete(); fHistAbso=0;
 255   fHistGcm2->Delete(); fHistGcm2=0;
 256   //
 257   if (fErrorCondition) DumpVolumes();
 258 }
 259
 260 //_______________________________________________________________________
 261 void AliLego::Copy(TObject&) const
 262 {
 263   //
 264   // Copy *this onto lego -- not implemented
 265   //
 266   Fatal("Copy","Not implemented!\n");
 267 }
 268
 269 //_______________________________________________________________________
 270 void AliLego::StepManager()
 271 {
 272   //
 273   // called from AliRun::Stepmanager from gustep.
 274   // Accumulate the 3 parameters step by step
 275   //
 276   static Float_t t;
 277   Float_t a,z,dens,radl,absl;
 278   Int_t i, id, copy;
 279   const char* vol;
 280   static Float_t vect[3], dir[3];
 281
 282   TString tmp1, tmp2;
 283   copy = 1;
 284   id  = gMC->CurrentVolID(copy);
 285   vol = gMC->VolName(id);
 286   Float_t step  = gMC->TrackStep();
 287
 288   TLorentzVector pos, mom;
 289   gMC->TrackPosition(pos);
 290   gMC->TrackMomentum(mom);
 291
 292   Int_t status = 0;
 293   if (gMC->IsTrackEntering()) status = 1;
 294   if (gMC->IsTrackExiting())  status = 2;
 295
 296   if (! fStepBack) {
 297     //      printf("\n volume %s %d", vol, status);
 298     //
 299     // Normal Forward stepping
 300     //
 301     if (fDebug) {
 302       //          printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
 303
 304       //
 305       // store volume if different from previous
 306       //
 307
 308           TClonesArray &lvols = *fVolumesFwd;
 309           if (fStepsForward > 0) {
 310         AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
 311         if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
 312                   fStepsForward -= 2;
 313                   fVolumesFwd->RemoveAt(fStepsForward);
 314                   fVolumesFwd->RemoveAt(fStepsForward+1);
 315         }
 316           }
 317
 318           new(lvols[fStepsForward++])
 319         AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
 320
 321     } // Debug
 322     //
 323     // Get current material properties
 324
 325     gMC->CurrentMaterial(a,z,dens,radl,absl);
 326
 327     if (z < 1) return;
 328
 329     // --- See if we have to stop now
 330     if (TMath::Abs(pos[2]) > fGener->ZMax()  ||
 331         pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
 332       if (!gMC->IsNewTrack()) {
 333         // Not the first step, add past contribution
 334         fTotAbso += t/absl;
 335         fTotRadl += t/radl;
 336         fTotGcm2 += t*dens;
 337         if (fDebug) {
 338           //
 339           //  generate "mirror" particle flying back
 340           //
 341           fStepsBackward = fStepsForward;
 342
 343           Float_t pmom[3], orig[3];
 344           Float_t polar[3] = {0.,0.,0.};
 345           Int_t ntr;
 346           pmom[0] = -dir[0];
 347           pmom[1] = -dir[1];
 348           pmom[2] = -dir[2];
 349           orig[0] =  vect[0];
 350           orig[1] =  vect[1];
 351           orig[2] =  vect[2];
 352
 353           gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(),
 354                            0, pmom, orig, polar, 0., kPNoProcess, ntr);
 355         } // debug
 356
 357       } // not a new track !
 358
 359       if (fDebug) fStepBack = 1;
 360       gMC->StopTrack();
 361       return;
 362     } // outside scoring region ?
 363
 364     // --- See how long we have to go
 365     for(i=0;i<3;++i) {
 366       vect[i]=pos[i];
 367       dir[i]=mom[i];
 368     }
 369
 370     t  = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
 371
 372     if(step) {
 373
 374       fTotAbso += step/absl;
 375       fTotRadl += step/radl;
 376       fTotGcm2 += step*dens;
 377     }
 378
 379   } else {
 380       if (fDebug) {
 381         //
 382         // Geometry debugging
 383         // Fly back and compare volume sequence
 384         //
 385         TClonesArray &lvols = *fVolumesBwd;
 386         if (fStepsBackward < fStepsForward) {
 387               AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
 388               if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
 389             fStepsBackward += 2;
 390             fVolumesBwd->RemoveAt(fStepsBackward-1);
 391             fVolumesBwd->RemoveAt(fStepsBackward-2);
 392               }
 393         }
 394
 395         fStepsBackward--;
 396         //        printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
 397         if (fStepsBackward < 0) {
 398               gMC->StopTrack();
 399               fStepBack = 0;
 400               return;
 401         }
 402
 403         new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
 404
 405         AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
 406         if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
 407           {
 408             printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
 409                    fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
 410             fErrorCondition = 1;
 411           }
 412       } // Debug
 413   } // bwd/fwd
 414 }
 415
 416 //_______________________________________________________________________
 417 void AliLego::DumpVolumes()
 418 {
 419   //
 420   // Dump volume sequence in case of error
 421   //
 422   printf("\n Dumping Volume Sequence:");
 423   printf("\n ==============================================");
 424
 425   for (Int_t i = fStepsForward-1; i>=0; i--)
 426     {
 427       AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
 428       AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
 429       if (tmp1)
 430         printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
 431                , i,
 432                tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
 433                tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
 434       if (tmp2 && i>= fStepsBackward)
 435         printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
 436                , i,
 437                tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
 438                tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
 439
 440       printf("\n ............................................................................\n");
 441     }
 442   printf("\n ==============================================\n");
 443 }