PHOS/AliPHOSvFast.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 // Implementation of the PHOS manager class for fast simulations
  20 // Tracks particles until the reach a grossly designed PHOS module
  21 // Modify the particles property (momentum, energy, type) according to
  22 //  the PHOS response function. The result is called a virtual reconstructed
  23 //  particle.
  24 //
  25 //*-- Author: Yves Schutz (SUBATECH)
  26
  27 // --- ROOT system ---
  28
  29 #include "TBRIK.h"
  30 #include "TNode.h"
  31 #include "TParticle.h"
  32 #include "TGeometry.h"
  33
  34 // --- Standard library ---
  35
  36 // --- AliRoot header files ---
  37 #include "AliPHOSFastRecParticle.h"
  38 #include "AliPHOSvFast.h"
  39 #include "AliPHOSGetter.h"
  40 #include "AliRun.h"
  41
  42 ClassImp(AliPHOSvFast)
  43
  44 //____________________________________________________________________________
  45 AliPHOSvFast::AliPHOSvFast() : AliPHOS()
  46 {
  47   // default ctor : initialize data member
  48    fBigBoxX = 0. ;
  49    fBigBoxY = 0. ;
  50    fBigBoxZ = 0. ;
  51    fFastRecParticles = 0 ;
  52    fNRecParticles = 0 ;
  53    fRan = 0 ;
  54    fResPara1 = 0. ;
  55    fResPara2 = 0. ;
  56    fResPara3 = 0. ;
  57    fPosParaA0 = 0. ;
  58    fPosParaA1 = 0. ;
  59    fPosParaB0 = 0. ;
  60    fPosParaB1 = 0. ;
  61    fPosParaB2 = 0. ;
  62
  63 }
  64
  65 //____________________________________________________________________________
  66 AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
  67   AliPHOS(name,title)
  68 {
  69   // ctor
  70
  71
  72   // create the Getter
  73   AliPHOSGetter::GetInstance(gDirectory->GetName(), 0) ;
  74
  75   SetBigBox(0, GetGeometry()->GetOuterBoxSize(0) ) ;
  76   SetBigBox(1, GetGeometry()->GetOuterBoxSize(3) + GetGeometry()->GetCPVBoxSize(1) ) ;
  77   SetBigBox(2, GetGeometry()->GetOuterBoxSize(2) );
  78
  79   fNRecParticles = 0 ;
  80   fFastRecParticles = new AliPHOSFastRecParticle::FastRecParticlesList("AliPHOSFastRecParticle", 100) ;
  81
  82   fResPara1 = 0.030 ;    // GeV
  83   fResPara2 = 0.00003 ;
  84   fResPara3 = 0.00001 ;
  85
  86   fPosParaA0 = 2.87 ;    // mm
  87   fPosParaA1 = -0.0975 ;
  88   fPosParaB0 = 0.257 ;
  89   fPosParaB1 = 0.137 ;
  90   fPosParaB2 = 0.00619 ;
  91 }
  92
  93 //____________________________________________________________________________
  94 AliPHOSvFast::AliPHOSvFast(const AliPHOSvFast & pv):
  95   AliPHOS(pv.GetName(), pv.GetTitle())
  96 {
  97   // cpy ctor
  98
  99
 100   // create the Getter
 101   AliPHOSGetter::GetInstance(gDirectory->GetName(), 0) ;
 102
 103   SetBigBox(0, pv.fBigBoxX) ;
 104   SetBigBox(1, pv.fBigBoxY) ;
 105   SetBigBox(2, pv.fBigBoxZ) ;
 106
 107   fNRecParticles = pv.fNRecParticles;
 108   fFastRecParticles = new TClonesArray( *(pv.fFastRecParticles) ) ;
 109
 110   fResPara1 = pv.fResPara1 ;
 111   fResPara2 = pv.fResPara2 ;
 112   fResPara3 = pv.fResPara3 ;
 113
 114   fPosParaA0 = pv.fPosParaA0 ;
 115   fPosParaA1 = pv.fPosParaA1 ;
 116   fPosParaB0 = pv.fPosParaB0 ;
 117   fPosParaB1 = pv.fPosParaB1 ;
 118   fPosParaB2 = pv.fPosParaB2 ;
 119 }
 120
 121 //____________________________________________________________________________
 122 AliPHOSvFast::~AliPHOSvFast()
 123 {
 124   // dtor
 125
 126   fFastRecParticles->Delete() ;
 127   delete fFastRecParticles ;
 128   fFastRecParticles = 0 ;
 129
 130 }
 131
 132 //____________________________________________________________________________
 133 void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
 134 {
 135   // Add a virtually reconstructed particle to the list
 136
 137   new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
 138   fNRecParticles++ ;
 139 }
 140
 141 //____________________________________________________________________________
 142 void AliPHOSvFast::BuildGeometry()
 143 {
 144   // Build the PHOS geometry for the ROOT display
 145    //BEGIN_HTML
 146   /*
 147     <H2>
 148      PHOS FAST in ALICE displayed by root
 149     </H2>
 150     <H4> All Views </H4>
 151     <P>
 152     <CENTER>
 153     <IMG Align=BOTTOM ALT="Fast All Views" SRC="../images/AliPHOSvFastAllViews.gif">
 154     </CENTER></P>
 155     <H4> Front View </H4>
 156     <P>
 157     <CENTER>
 158     <IMG Align=BOTTOM ALT="Fast Front View" SRC="../images/AliPHOSvFastFrontView.gif">
 159     </CENTER></P>
 160   */
 161   //END_HTML
 162
 163   const Int_t kColorPHOS = kRed ;
 164
 165   Double_t const kRADDEG = 180.0 / kPI ;
 166
 167   new TBRIK( "BigBox", "PHOS box", "void", GetBigBox(0)/2,
 168              GetBigBox(1)/2,
 169              GetBigBox(2)/2 );
 170
 171   // position PHOS into ALICE
 172
 173   Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
 174   Int_t number = 988 ;
 175   Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * GetGeometry()->GetIPtoCrystalSurface() ) ) ;
 176   pphi *= kRADDEG ;
 177   TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
 178
 179   char * nodename = new char[20] ;
 180   char * rotname  = new char[20] ;
 181
 182   for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
 183    Float_t angle = pphi * 2 * ( i - GetGeometry()->GetNModules() / 2.0 - 0.5 ) ;
 184    sprintf(rotname, "%s%d", "rot", number++) ;
 185    new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
 186    top->cd();
 187    sprintf(nodename,"%s%d", "Module", i) ;
 188    Float_t x =  r * TMath::Sin( angle / kRADDEG ) ;
 189    Float_t y = -r * TMath::Cos( angle / kRADDEG ) ;
 190    TNode * bigboxnode = new TNode(nodename, nodename, "BigBox", x, y, 0, rotname ) ;
 191    bigboxnode->SetLineColor(kColorPHOS) ;
 192    fNodes->Add(bigboxnode) ;
 193   }
 194   delete[] nodename ;
 195   delete[] rotname ;
 196 }
 197
 198 //____________________________________________________________________________
 199 void AliPHOSvFast::CreateGeometry()
 200 {
 201   // Create the geometry for GEANT
 202
 203   AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
 204
 205   if ( phostmp == NULL ) {
 206
 207     fprintf(stderr, "PHOS detector not found!\n") ;
 208     return ;
 209
 210   }
 211
 212   // Get pointer to the array containing media indeces
 213   Int_t *idtmed = fIdtmed->GetArray() - 699 ;
 214
 215   Float_t bigbox[3] ;
 216   bigbox[0] =   GetBigBox(0) / 2.0 ;
 217   bigbox[1] =   GetBigBox(1) / 2.0 ;
 218   bigbox[2] =   GetBigBox(2) / 2.0 ;
 219
 220   gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
 221
 222   // --- Position  PHOS mdules in ALICE setup ---
 223
 224   Int_t idrotm[99] ;
 225   Double_t const kRADDEG = 180.0 / kPI ;
 226
 227   for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
 228
 229     Float_t angle = GetGeometry()->GetPHOSAngle(i) ;
 230     AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
 231
 232     Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
 233
 234     Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
 235     Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
 236     gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
 237
 238   } // for GetNModules
 239
 240 }
 241
 242
 243 //____________________________________________________________________________
 244 void AliPHOSvFast::Init(void)
 245 {
 246   // Prints out an information message
 247
 248   Int_t i;
 249
 250   printf("\n");
 251   for(i=0;i<35;i++) printf("*");
 252   printf(" FAST PHOS_INIT ");
 253   for(i=0;i<35;i++) printf("*");
 254   printf("\n");
 255
 256   // Here the PHOS initialisation code (if any!)
 257
 258   for(i=0;i<80;i++) printf("*");
 259   printf("\n");
 260
 261 }
 262
 263 //___________________________________________________________________________
 264 Float_t AliPHOSvFast::GetBigBox(Int_t index) const
 265 {
 266   // Get the X, Y or Z dimension of the box describing a PHOS module
 267
 268   Float_t rv = 0 ;
 269
 270   switch (index) {
 271   case 0:
 272     rv = fBigBoxX ;
 273     break ;
 274   case 1:
 275      rv = fBigBoxY ;
 276     break ;
 277   case 2:
 278      rv = fBigBoxZ ;
 279     break ;
 280  }
 281   return rv ;
 282 }
 283
 284 //___________________________________________________________________________
 285 void AliPHOSvFast::MakeBranch(Option_t* opt, const char *file)
 286 {
 287   // Create new branch in the current reconstructed Root Tree
 288
 289   AliDetector::MakeBranch(opt,file) ;
 290
 291   char branchname[10];
 292   sprintf(branchname,"%s",GetName());
 293   const char *cd = strstr(opt,"R");
 294
 295   if (fFastRecParticles && gAlice->TreeR() && cd) {
 296     MakeBranchInTree(gAlice->TreeR(),
 297                      branchname, &fFastRecParticles, fBufferSize, file);
 298   }
 299 }
 300
 301 //____________________________________________________________________________
 302 Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
 303 {
 304   // Smears the energy according to the energy dependent energy resolution.
 305   // A gaussian distribution is assumed
 306
 307   Double_t sigma  = SigmaE(energy) ;
 308   return  fRan.Gaus(energy, sigma) ;
 309 }
 310
 311 //____________________________________________________________________________
 312 TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi)
 313 {
 314   // Smears the impact position according to the energy dependent position resolution
 315   // A gaussian position distribution is assumed
 316
 317   TVector3 newpos ;
 318   Double_t sigma = SigmaP( energy, theta*180./TMath::Pi() ) ;
 319   Double_t x = fRan.Gaus( pos.X(), sigma ) ;
 320   sigma = SigmaP( energy, phi*180./TMath::Pi() ) ;
 321   Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
 322   Double_t y = pos.Y() ;
 323
 324   newpos.SetX(x) ;
 325   newpos.SetY(y) ;
 326   newpos.SetZ(z) ;
 327
 328   return newpos ;
 329 }
 330
 331 //____________________________________________________________________________
 332 void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
 333 {
 334   // Modify the primary particle properties according
 335   //  1. the response function of PHOS
 336   //  2. the performance of the EMC+PPSD setup
 337
 338   Int_t type = MakeType( rp ) ;
 339   rp.SetType(type) ;
 340
 341
 342   // get the detected energy
 343
 344   TLorentzVector momentum ;
 345   rp.Momentum(momentum) ;
 346   Double_t kineticenergy = TMath::Sqrt( TMath::Power(momentum.E(), 2) - TMath::Power(rp.GetMass(), 2) ) ;
 347   Double_t modifiedkineticenergy = MakeEnergy(kineticenergy ) ;
 348   Double_t modifiedenergy = TMath::Sqrt( TMath::Power(modifiedkineticenergy, 2)
 349                                          + TMath::Power( rp.GetMass(), 2) ) ;
 350
 351   // get the angle of incidence
 352
 353   Double_t incidencetheta = 90. * TMath::Pi() /180 - rp.Theta() ;
 354   Double_t incidencephi   = ( 270 + GetGeometry()->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;
 355
 356   // get the detected direction
 357
 358   TVector3 modifiedposition = MakePosition(kineticenergy, pos, incidencetheta, incidencephi) ;
 359   modifiedposition *= modifiedkineticenergy / modifiedposition.Mag() ;
 360
 361   // Set the modified 4-momentum of the reconstructed particle
 362
 363   rp.SetMomentum(modifiedposition.X(), modifiedposition.Y(), modifiedposition.Z(), modifiedenergy) ;
 364
 365  }
 366
 367 //____________________________________________________________________________
 368 Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
 369 {
 370   // Generate a particle type using the performance of the EMC+PPSD setup
 371
 372   Int_t rv =   AliPHOSFastRecParticle::kUNDEFINED ;
 373   Int_t charge = (Int_t)rp.GetPDG()->Charge() ;
 374   Int_t test ;
 375   Float_t ran ;
 376   if ( charge != 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) )
 377     test = - 1 ;
 378   else
 379     test = rp.GetPdgCode() ;
 380
 381   Fatal("MakeType", "SHOULD NOT BE USED until values of probabilities are properly set ") ;
 382
 383   switch (test) {
 384
 385   case 22:    // it's a photon              // NB: ALL VALUES SHOLD BE CHECKED !!!!
 386     ran = fRan.Rndm() ;
 387     if( ran <= 0.9498 )
 388       rv =  AliPHOSFastRecParticle::kNEUTRALHAFAST ;
 389     else
 390       rv =  AliPHOSFastRecParticle::kNEUTRALEMFAST ;
 391     break ;
 392
 393   case 2112:  // it's a neutron
 394     ran = fRan.Rndm() ;
 395     if ( ran <= 0.9998 )
 396       rv =  AliPHOSFastRecParticle::kNEUTRALHASLOW ;
 397     else
 398       rv = AliPHOSFastRecParticle::kNEUTRALEMSLOW ;
 399     break ;
 400
 401   case -2112: // it's a anti-neutron
 402     ran = fRan.Rndm() ;
 403     if ( ran <= 0.9984 )
 404       rv =  AliPHOSFastRecParticle::kNEUTRALHASLOW ;
 405     else
 406       rv =  AliPHOSFastRecParticle::kNEUTRALEMSLOW ;
 407     break ;
 408
 409   case 11:    // it's a electron
 410     ran = fRan.Rndm() ;
 411     if ( ran <= 0.9996 )
 412       rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ;
 413     else
 414       rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ;
 415     break;
 416
 417   case -11:   // it's a positon
 418     ran = fRan.Rndm() ;
 419     if ( ran <= 0.9996 )
 420       rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ;
 421     else
 422       rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ;
 423     break;
 424
 425   case -1:    // it's a charged
 426     ran = fRan.Rndm() ;
 427     if ( ran <= 0.9996 )
 428       rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ;
 429     else
 430       rv =  AliPHOSFastRecParticle::kNEUTRALHAFAST ;
 431
 432     break ;
 433   }
 434
 435
 436   return rv ;
 437 }
 438
 439 //___________________________________________________________________________
 440 void AliPHOSvFast::ResetPoints()
 441 {
 442   // This overloads the method in AliDetector
 443
 444   ResetFastRecParticles() ;
 445 }
 446
 447 //___________________________________________________________________________
 448 void AliPHOSvFast::ResetFastRecParticles()
 449 {
 450   // Resets the list of virtual reconstructed particles
 451
 452   if (fFastRecParticles)
 453     fFastRecParticles->Clear() ;
 454   fNRecParticles = 0 ;
 455 }
 456
 457 //___________________________________________________________________________
 458 void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
 459 {
 460   // Set the size of the Box describing a PHOS module
 461
 462   switch (index) {
 463   case 0:
 464     fBigBoxX = value ;
 465     break ;
 466   case 1:
 467     fBigBoxY = value ;
 468     break ;
 469   case 2:
 470     fBigBoxZ = value ;
 471     break ;
 472  }
 473
 474 }
 475
 476 //____________________________________________________________________________
 477 Double_t AliPHOSvFast::SigmaE(Double_t energy)
 478 {
 479   // Calculates the energy dependent energy resolution
 480
 481   Double_t rv = -1 ;
 482
 483   rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
 484                + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
 485                + TMath::Power(fResPara3, 2) ) ;
 486
 487   return rv * energy ;
 488 }
 489
 490 //____________________________________________________________________________
 491 Double_t AliPHOSvFast::SigmaP(Double_t energy, Double_t incidence)
 492 {
 493   // Calculates the energy dependent position resolution
 494
 495   Double_t paraA = fPosParaA0 + fPosParaA1 * incidence ;
 496   Double_t paraB = fPosParaB0 + fPosParaB1 * incidence + fPosParaB2 * incidence * incidence ;
 497
 498   return ( paraA / TMath::Sqrt(energy) + paraB ) * 0.1   ; // in cm
 499 }
 500
 501 //____________________________________________________________________________
 502 void AliPHOSvFast::StepManager(void)
 503 {
 504   // Only verifies if the particle reaches PHOS and stops the tracking
 505
 506   TLorentzVector lv ;
 507   gMC->TrackPosition(lv) ;
 508   TVector3 pos = lv.Vect() ;
 509   Int_t modid  ;
 510   gMC->CurrentVolID(modid);
 511
 512   Float_t energy = gMC->Etot() ; //Total energy of current track
 513
 514   //Calculating mass of current particle
 515   TDatabasePDG * pdg = TDatabasePDG::Instance() ;
 516   TParticlePDG * partPDG = pdg->GetParticle(gMC->TrackPid()) ;
 517   Float_t mass = partPDG->Mass() ;
 518
 519   if(energy > mass){
 520     pos.SetMag(TMath::Sqrt(energy*energy-mass*mass)) ;
 521     TLorentzVector pTrack(pos, energy) ;
 522
 523     TParticle * part = new TParticle(gMC->TrackPid(), 0,-1,-1,-1,-1, pTrack, lv)  ;
 524
 525     AliPHOSFastRecParticle rp(*part) ;
 526
 527     // Adds the response of PHOS to the particle
 528     MakeRecParticle(modid, pos, rp) ;
 529
 530     // add the `track' particle to the FastRecParticles list
 531
 532     AddRecParticle(rp) ;
 533
 534     part->Delete() ;
 535   }
 536   // stop the track as soon PHOS is reached
 537
 538   gMC->StopTrack() ;
 539
 540 }
 541