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