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