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