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 {
  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          = AliDebugLevel();
 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          = AliDebugLevel();
 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) ToAliDebug(1, 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) ToAliError(DumpVolumes());
 259 }
 260
 261 //_______________________________________________________________________
 262 void AliLego::Copy(TObject&) const
 263 {
 264   //
 265   // Copy *this onto lego -- not implemented
 266   //
 267   AliFatal("Not implemented!");
 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                 if (absl) fTotAbso += t/absl;
 336                 if (radl) 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             if (absl) fTotAbso += step/absl;
 376             if (radl) 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                 AliWarning(Form("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 }