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