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