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