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