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