Effective C++ corrections (T.Pocheptsov)
[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$
43fbaae1 21 * Revision 1.29 2005/05/28 14:19:05 schutz
22 * Compilation warnings fixed by T.P.
23 *
702ab87e 24 */
25
3decf5cb 26//_________________________________________________________________________
b2a60966 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)
3decf5cb 34
35// --- ROOT system ---
3a4295ee 36
88cb7938 37#include <TBRIK.h>
88cb7938 38#include <TGeometry.h>
39#include <TNode.h>
40#include <TParticle.h>
e957fea8 41#include "TClonesArray.h"
88cb7938 42#include <TVirtualMC.h>
3decf5cb 43
44// --- Standard library ---
45
3decf5cb 46// --- AliRoot header files ---
3a4295ee 47#include "AliPHOSFastRecParticle.h"
88cb7938 48#include "AliPHOSGeometry.h"
49#include "AliPHOSLoader.h"
3decf5cb 50#include "AliPHOSvFast.h"
3decf5cb 51#include "AliRun.h"
3decf5cb 52
53ClassImp(AliPHOSvFast)
54
43fbaae1 55AliPHOSvFast::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.)
3decf5cb 70{
3a4295ee 71 // default ctor : initialize data member
3decf5cb 72}
73
74//____________________________________________________________________________
75AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
43fbaae1 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)
3decf5cb 91{
b2a60966 92 // ctor
88cb7938 93 // create the Loader
7a9d98f9 94 SetBigBox(0, GetGeometry()->GetOuterBoxSize(0) ) ;
95 SetBigBox(1, GetGeometry()->GetOuterBoxSize(3) + GetGeometry()->GetCPVBoxSize(1) ) ;
96 SetBigBox(2, GetGeometry()->GetOuterBoxSize(2) );
3decf5cb 97}
98
99//____________________________________________________________________________
43fbaae1 100AliPHOSvFast::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//____________________________________________________________________________
3decf5cb 123AliPHOSvFast::~AliPHOSvFast()
124{
b2a60966 125 // dtor
3decf5cb 126
a73f33f0 127 fFastRecParticles->Delete() ;
128 delete fFastRecParticles ;
129 fFastRecParticles = 0 ;
3decf5cb 130
131}
132
133//____________________________________________________________________________
702ab87e 134void AliPHOSvFast::Copy(TObject & base)const
780fda6d 135{
702ab87e 136 TObject::Copy(base) ;
137 AliPHOS::Copy(base) ;
138 AliPHOSvFast &fast = static_cast<AliPHOSvFast &>(base);
780fda6d 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//____________________________________________________________________________
a73f33f0 158void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
159{
b2a60966 160 // Add a virtually reconstructed particle to the list
161
f96d3dc7 162 new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
163 fNRecParticles++ ;
3decf5cb 164}
165
166//____________________________________________________________________________
167void AliPHOSvFast::BuildGeometry()
168{
3decf5cb 169 // Build the PHOS geometry for the ROOT display
b2a60966 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
3decf5cb 188 const Int_t kColorPHOS = kRed ;
189
a8c47ab6 190 Double_t const kRADDEG = 180.0 / TMath::Pi() ;
3decf5cb 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
fa7cce36 198 Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
3decf5cb 199 Int_t number = 988 ;
fa7cce36 200 Float_t pphi = TMath::ATan( GetBigBox(0) / ( 2.0 * GetGeometry()->GetIPtoCrystalSurface() ) ) ;
3decf5cb 201 pphi *= kRADDEG ;
202 TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
203
204 char * nodename = new char[20] ;
205 char * rotname = new char[20] ;
206
fa7cce36 207 for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
208 Float_t angle = pphi * 2 * ( i - GetGeometry()->GetNModules() / 2.0 - 0.5 ) ;
3decf5cb 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//____________________________________________________________________________
224void AliPHOSvFast::CreateGeometry()
225{
b2a60966 226 // Create the geometry for GEANT
227
3decf5cb 228 AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
b2a60966 229
3decf5cb 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] ;
a8c47ab6 250 Double_t const kRADDEG = 180.0 / TMath::Pi() ;
3decf5cb 251
fa7cce36 252 for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
3decf5cb 253
fa7cce36 254 Float_t angle = GetGeometry()->GetPHOSAngle(i) ;
3decf5cb 255 AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
256
fa7cce36 257 Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
3decf5cb 258
259 Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
260 Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
3decf5cb 261 gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
262
263 } // for GetNModules
264
265}
266
267
268//____________________________________________________________________________
269void AliPHOSvFast::Init(void)
270{
b2a60966 271 // Prints out an information message
272
3decf5cb 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//___________________________________________________________________________
3a4295ee 289Float_t AliPHOSvFast::GetBigBox(Int_t index) const
3decf5cb 290{
b2a60966 291 // Get the X, Y or Z dimension of the box describing a PHOS module
292
3decf5cb 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}
3decf5cb 308//___________________________________________________________________________
88cb7938 309
310void AliPHOSvFast::MakeBranch(Option_t* opt)
3decf5cb 311{
b2a60966 312 // Create new branch in the current reconstructed Root Tree
88cb7938 313 AliDetector::MakeBranch(opt);
3a4295ee 314 const char *cd = strstr(opt,"R");
3decf5cb 315
88cb7938 316 if (fFastRecParticles && fLoader->TreeR() && cd) {
317 MakeBranchInTree(fLoader->TreeR(), GetName(), &fFastRecParticles, fBufferSize, 0);
3decf5cb 318 }
319}
a73f33f0 320//____________________________________________________________________________
88cb7938 321
fc7e2f43 322Double_t AliPHOSvFast::MakeEnergy(Double_t energy)
a73f33f0 323{
b2a60966 324 // Smears the energy according to the energy dependent energy resolution.
325 // A gaussian distribution is assumed
326
b1aa729d 327 Double_t sigma = SigmaE(energy) ;
328 return fRan.Gaus(energy, sigma) ;
a73f33f0 329}
a73f33f0 330//____________________________________________________________________________
88cb7938 331
fc7e2f43 332TVector3 AliPHOSvFast::MakePosition(Double_t energy, TVector3 pos, Double_t theta, Double_t phi)
a73f33f0 333{
b2a60966 334 // Smears the impact position according to the energy dependent position resolution
335 // A gaussian position distribution is assumed
336
b1aa729d 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}
a73f33f0 350
b1aa729d 351//____________________________________________________________________________
fc7e2f43 352void AliPHOSvFast::MakeRecParticle(Int_t modid, TVector3 pos, AliPHOSFastRecParticle & rp)
b1aa729d 353{
b2a60966 354 // Modify the primary particle properties according
355 // 1. the response function of PHOS
356 // 2. the performance of the EMC+PPSD setup
357
b1aa729d 358 Int_t type = MakeType( rp ) ;
a73f33f0 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 ) ;
b1aa729d 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() ;
fa7cce36 374 Double_t incidencephi = ( 270 + GetGeometry()->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;
b1aa729d 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
a73f33f0 385 }
386
387//____________________________________________________________________________
b1aa729d 388Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
a73f33f0 389{
b2a60966 390 // Generate a particle type using the performance of the EMC+PPSD setup
391
9ec91567 392 Int_t rv = AliPHOSFastRecParticle::kUNDEFINED ;
b1aa729d 393 Int_t charge = (Int_t)rp.GetPDG()->Charge() ;
394 Int_t test ;
395 Float_t ran ;
f96d3dc7 396 if ( charge != 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) )
b1aa729d 397 test = - 1 ;
398 else
399 test = rp.GetPdgCode() ;
400
a4765fe4 401 Fatal("MakeType", "SHOULD NOT BE USED until values of probabilities are properly set ") ;
402 // NB: ALL VALUES SHOULD BE CHECKED !!!!
b1aa729d 403 switch (test) {
404
3a4295ee 405 case 22: // it's a photon // NB: ALL VALUES SHOLD BE CHECKED !!!!
b1aa729d 406 ran = fRan.Rndm() ;
3a4295ee 407 if( ran <= 0.9498 )
408 rv = AliPHOSFastRecParticle::kNEUTRALHAFAST ;
409 else
410 rv = AliPHOSFastRecParticle::kNEUTRALEMFAST ;
b1aa729d 411 break ;
412
413 case 2112: // it's a neutron
414 ran = fRan.Rndm() ;
415 if ( ran <= 0.9998 )
3a4295ee 416 rv = AliPHOSFastRecParticle::kNEUTRALHASLOW ;
b1aa729d 417 else
3a4295ee 418 rv = AliPHOSFastRecParticle::kNEUTRALEMSLOW ;
b1aa729d 419 break ;
3a4295ee 420
b1aa729d 421 case -2112: // it's a anti-neutron
422 ran = fRan.Rndm() ;
423 if ( ran <= 0.9984 )
3a4295ee 424 rv = AliPHOSFastRecParticle::kNEUTRALHASLOW ;
b1aa729d 425 else
3a4295ee 426 rv = AliPHOSFastRecParticle::kNEUTRALEMSLOW ;
b1aa729d 427 break ;
428
429 case 11: // it's a electron
430 ran = fRan.Rndm() ;
431 if ( ran <= 0.9996 )
3a4295ee 432 rv = AliPHOSFastRecParticle::kCHARGEDEMFAST ;
b1aa729d 433 else
3a4295ee 434 rv = AliPHOSFastRecParticle::kCHARGEDHAFAST ;
b1aa729d 435 break;
436
f96d3dc7 437 case -11: // it's a positon
b1aa729d 438 ran = fRan.Rndm() ;
439 if ( ran <= 0.9996 )
3a4295ee 440 rv = AliPHOSFastRecParticle::kCHARGEDEMFAST ;
b1aa729d 441 else
3a4295ee 442 rv = AliPHOSFastRecParticle::kCHARGEDHAFAST ;
b1aa729d 443 break;
444
445 case -1: // it's a charged
446 ran = fRan.Rndm() ;
447 if ( ran <= 0.9996 )
3a4295ee 448 rv = AliPHOSFastRecParticle::kCHARGEDHAFAST ;
b1aa729d 449 else
3a4295ee 450 rv = AliPHOSFastRecParticle::kNEUTRALHAFAST ;
b1aa729d 451
452 break ;
453 }
454
455
a73f33f0 456 return rv ;
457}
458
3decf5cb 459//___________________________________________________________________________
f96d3dc7 460void AliPHOSvFast::ResetPoints()
461{
b2a60966 462 // This overloads the method in AliDetector
463
f96d3dc7 464 ResetFastRecParticles() ;
465}
466
467//___________________________________________________________________________
468void AliPHOSvFast::ResetFastRecParticles()
469{
b2a60966 470 // Resets the list of virtual reconstructed particles
471
f96d3dc7 472 if (fFastRecParticles)
473 fFastRecParticles->Clear() ;
474 fNRecParticles = 0 ;
475}
476
477//___________________________________________________________________________
3decf5cb 478void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
479{
b2a60966 480 // Set the size of the Box describing a PHOS module
481
3decf5cb 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//____________________________________________________________________________
a73f33f0 497Double_t AliPHOSvFast::SigmaE(Double_t energy)
498{
b2a60966 499 // Calculates the energy dependent energy resolution
500
a73f33f0 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//____________________________________________________________________________
7a9d98f9 511Double_t AliPHOSvFast::SigmaP(Double_t energy, Double_t incidence)
b1aa729d 512{
b2a60966 513 // Calculates the energy dependent position resolution
514
b1aa729d 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//____________________________________________________________________________
3decf5cb 522void AliPHOSvFast::StepManager(void)
523{
b2a60966 524 // Only verifies if the particle reaches PHOS and stops the tracking
3a4295ee 525
a73f33f0 526 TLorentzVector lv ;
527 gMC->TrackPosition(lv) ;
b1aa729d 528 TVector3 pos = lv.Vect() ;
529 Int_t modid ;
530 gMC->CurrentVolID(modid);
a73f33f0 531
3a4295ee 532 Float_t energy = gMC->Etot() ; //Total energy of current track
3decf5cb 533
3a4295ee 534 //Calculating mass of current particle
535 TDatabasePDG * pdg = TDatabasePDG::Instance() ;
536 TParticlePDG * partPDG = pdg->GetParticle(gMC->TrackPid()) ;
537 Float_t mass = partPDG->Mass() ;
a73f33f0 538
3a4295ee 539 if(energy > mass){
540 pos.SetMag(TMath::Sqrt(energy*energy-mass*mass)) ;
541 TLorentzVector pTrack(pos, energy) ;
a73f33f0 542
3a4295ee 543 TParticle * part = new TParticle(gMC->TrackPid(), 0,-1,-1,-1,-1, pTrack, lv) ;
544
545 AliPHOSFastRecParticle rp(*part) ;
b1aa729d 546
3a4295ee 547 // Adds the response of PHOS to the particle
548 MakeRecParticle(modid, pos, rp) ;
549
550 // add the `track' particle to the FastRecParticles list
b1aa729d 551
3a4295ee 552 AddRecParticle(rp) ;
a73f33f0 553
3a4295ee 554 part->Delete() ;
555 }
a73f33f0 556 // stop the track as soon PHOS is reached
b1aa729d 557
3decf5cb 558 gMC->StopTrack() ;
559
560}
561