Fast simulation ready to be tested
[u/mrichter/AliRoot.git] / 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 //_________________________________________________________________________
17 // Manager class for PHOS version for fast simulations
18 //*-- Author : Y. Schutz SUBATECH 
19 //////////////////////////////////////////////////////////////////////////////
20
21 // --- ROOT system ---
22
23 #include "TBRIK.h"
24 #include "TNode.h"
25 #include "TParticle.h"
26
27 // --- Standard library ---
28
29 #include <cstdio>
30 #include <cassert>
31
32 // --- AliRoot header files ---
33
34 #include "AliPHOSvFast.h"
35 #include "AliRun.h"
36 #include "AliConst.h"
37
38 ClassImp(AliPHOSvFast)
39
40 //____________________________________________________________________________
41 AliPHOSvFast::AliPHOSvFast()
42 {
43   fFastRecParticles = 0 ; 
44   fNRecParticles = 0 ; 
45 }
46
47 //____________________________________________________________________________
48 AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
49   AliPHOS(name,title)
50 {
51   // gets an instance of the geometry parameters class  
52    
53   fGeom =  AliPHOSGeometry::GetInstance(title, "") ; 
54
55   if (fGeom->IsInitialized() ) 
56     cout << "AliPHOSvFast : PHOS geometry intialized for " << fGeom->GetName() << endl ;
57   else
58     cout << "AliPHOSvFast : PHOS geometry initialization failed !" << endl ;   
59   
60   SetBigBox(0, fGeom->GetOuterBoxSize(0) ) ;
61   SetBigBox(1, fGeom->GetOuterBoxSize(1) + fGeom->GetPPSDBoxSize(1) ) ; 
62   SetBigBox(2, fGeom->GetOuterBoxSize(0) ); 
63
64   fNRecParticles = 0 ; 
65   fFastRecParticles = new FastRecParticlesList("AliPHOSFastRecParticle", 100) ;
66
67   fResPara1 = 0.030 ;    // GeV
68   fResPara2 = 0.00003 ; 
69   fResPara3 = 0.00001 ; 
70
71   fPosParaA0 = 2.87 ;    // mm
72   fPosParaA1 = -0.0975 ;  
73   fPosParaB0 = 0.257 ;   
74   fPosParaB1 = 0.137 ; 
75   fPosParaB2 = 0.00619 ; 
76 }
77
78 //____________________________________________________________________________
79 AliPHOSvFast::~AliPHOSvFast()
80 {
81  
82   fFastRecParticles->Delete() ; 
83   delete fFastRecParticles ;
84   fFastRecParticles = 0 ; 
85
86 }
87
88 //____________________________________________________________________________
89 void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
90 {  
91   new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
92   fNRecParticles++ ; 
93 }
94
95 //____________________________________________________________________________
96 void AliPHOSvFast::BuildGeometry()
97 {
98
99   // Build the PHOS geometry for the ROOT display
100   
101   const Int_t kColorPHOS = kRed ;
102   
103   Double_t const kRADDEG = 180.0 / kPI ;
104   
105   new TBRIK( "BigBox", "PHOS box", "void", GetBigBox(0)/2, 
106              GetBigBox(1)/2, 
107              GetBigBox(2)/2 );
108   
109   // position PHOS into ALICE
110
111   Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
112   Int_t number = 988 ; 
113   Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * fGeom->GetIPtoOuterCoverDistance() ) ) ;
114   pphi *= kRADDEG ;
115   TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
116  
117   char * nodename = new char[20] ;  
118   char * rotname  = new char[20] ; 
119
120   for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { 
121    Float_t angle = pphi * 2 * ( i - fGeom->GetNModules() / 2.0 - 0.5 ) ;
122    sprintf(rotname, "%s%d", "rot", number++) ;
123    new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
124    top->cd();
125    sprintf(nodename,"%s%d", "Module", i) ;    
126    Float_t x =  r * TMath::Sin( angle / kRADDEG ) ;
127    Float_t y = -r * TMath::Cos( angle / kRADDEG ) ;
128    TNode * bigboxnode = new TNode(nodename, nodename, "BigBox", x, y, 0, rotname ) ;
129    bigboxnode->SetLineColor(kColorPHOS) ;
130    fNodes->Add(bigboxnode) ;
131   }
132   delete[] nodename ; 
133   delete[] rotname ; 
134 }
135
136 //____________________________________________________________________________
137 void AliPHOSvFast::CreateGeometry()
138 {
139
140   AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
141
142   if ( phostmp == NULL ) {
143     
144     fprintf(stderr, "PHOS detector not found!\n") ;
145     return ;
146     
147   }
148
149   // Get pointer to the array containing media indeces
150   Int_t *idtmed = fIdtmed->GetArray() - 699 ;
151   
152   Float_t bigbox[3] ; 
153   bigbox[0] =   GetBigBox(0) / 2.0 ;
154   bigbox[1] =   GetBigBox(1) / 2.0 ;
155   bigbox[2] =   GetBigBox(2) / 2.0 ;
156   
157   gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
158   
159   // --- Position  PHOS mdules in ALICE setup ---
160   
161   Int_t idrotm[99] ;
162   Double_t const kRADDEG = 180.0 / kPI ;
163   
164   for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
165     
166     Float_t angle = fGeom->GetPHOSAngle(i) ;
167     AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
168  
169     Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
170
171     Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
172     Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
173     gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
174  
175   } // for GetNModules
176
177 }
178
179
180 //____________________________________________________________________________
181 void AliPHOSvFast::Init(void)
182 {
183  
184   Int_t i;
185
186   printf("\n");
187   for(i=0;i<35;i++) printf("*");
188   printf(" FAST PHOS_INIT ");
189   for(i=0;i<35;i++) printf("*");
190   printf("\n");
191
192   // Here the PHOS initialisation code (if any!)
193
194   for(i=0;i<80;i++) printf("*");
195   printf("\n");
196   
197 }
198
199 //___________________________________________________________________________
200 Float_t AliPHOSvFast::GetBigBox(Int_t index)
201 {
202   Float_t rv = 0 ; 
203
204   switch (index) {
205   case 0:
206     rv = fBigBoxX ; 
207     break ; 
208   case 1:
209      rv = fBigBoxY ; 
210     break ; 
211   case 2:
212      rv = fBigBoxZ ; 
213     break ; 
214  }
215   return rv ; 
216 }
217
218 //___________________________________________________________________________
219 void AliPHOSvFast::MakeBranch(Option_t* opt)
220 {  
221   //
222   // Create a new branch in the current Root Tree
223  
224   AliDetector::MakeBranch(opt) ;
225   
226   char branchname[10];
227   sprintf(branchname,"%s",GetName());
228   char *cd = strstr(opt,"R");
229   
230   if (fFastRecParticles && gAlice->TreeR() && cd) {
231     gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
232   }
233 }
234
235 //____________________________________________________________________________
236 Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
237 {  
238   Double_t sigma  = SigmaE(energy) ; 
239   return  fRan.Gaus(energy, sigma) ;   
240 }
241
242 //____________________________________________________________________________
243 TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi)
244 {
245   TVector3 newpos ;
246   Double_t sigma = SigmaP( energy, theta*180./TMath::Pi() ) ;
247   Double_t x = fRan.Gaus( pos.X(), sigma ) ;
248   sigma = SigmaP( energy, phi*180./TMath::Pi() ) ;
249   Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
250   Double_t y = pos.Y() ; 
251   
252   newpos.SetX(x) ; 
253   newpos.SetY(y) ; 
254   newpos.SetZ(z) ; 
255               
256   return newpos ; 
257 }
258
259 //____________________________________________________________________________
260 void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
261 {
262   // get the detected type of particle
263  
264   Int_t type = MakeType( rp ) ;
265   rp.SetType(type) ;
266
267   
268   // get the detected energy
269
270   TLorentzVector momentum ;  
271   rp.Momentum(momentum) ; 
272   Double_t kineticenergy = TMath::Sqrt( TMath::Power(momentum.E(), 2) - TMath::Power(rp.GetMass(), 2) ) ; 
273   Double_t modifiedkineticenergy = MakeEnergy(kineticenergy ) ;
274   Double_t modifiedenergy = TMath::Sqrt( TMath::Power(modifiedkineticenergy, 2)  
275                                          + TMath::Power( rp.GetMass(), 2) ) ;
276  
277   // get the angle of incidence 
278   
279   Double_t incidencetheta = 90. * TMath::Pi() /180 - rp.Theta() ; 
280   Double_t incidencephi   = ( 270 + fGeom->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;   
281
282   // get the detected direction
283   
284   TVector3 modifiedposition = MakePosition(kineticenergy, pos, incidencetheta, incidencephi) ; 
285   modifiedposition *= modifiedkineticenergy / modifiedposition.Mag() ; 
286
287   // Set the modified 4-momentum of the reconstructed particle
288
289   rp.SetMomentum(modifiedposition.X(), modifiedposition.Y(), modifiedposition.Z(), modifiedenergy) ; 
290
291  }
292
293 //____________________________________________________________________________
294 Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
295 {
296   Int_t rv =  kUNDEFINED ;
297   Int_t charge = (Int_t)rp.GetPDG()->Charge() ;
298   Int_t test ; 
299   Float_t ran ; 
300   if ( charge != 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) ) 
301     test = - 1 ;
302   else
303     test = rp.GetPdgCode() ; 
304
305   switch (test) { 
306
307   case 22:    // it's a photon
308     ran = fRan.Rndm() ; 
309     if ( ran <= 0.5 )  // 50 % 
310       rv = kGAMMA ; 
311     else {
312       ran = fRan.Rndm() ; 
313       if( ran <= 0.9498 )
314         rv = kNEUTRALEM ; 
315       else
316         rv = kNEUTRALHADRON ; 
317     }     
318     break ; 
319
320   case 2112:  // it's a neutron
321     ran = fRan.Rndm() ; 
322     if ( ran <= 0.9998 )
323       rv = kNEUTRALHADRON ; 
324     else 
325       rv = kNEUTRALEM ; 
326     break ; 
327
328   case -2112: // it's a anti-neutron
329     ran = fRan.Rndm() ; 
330     if ( ran <= 0.9984 )
331       rv = kNEUTRALHADRON ; 
332     else 
333       rv = kNEUTRALEM ; 
334     break ; 
335     
336   case 11:    // it's a electron
337     ran = fRan.Rndm() ; 
338     if ( ran <= 0.9996 )
339       rv = kELECTRON ; 
340     else 
341       rv = kCHARGEDHADRON ; 
342     break; 
343
344   case -11:   // it's a positon
345     ran = fRan.Rndm() ; 
346     if ( ran <= 0.9996 )
347       rv = kELECTRON ; 
348     else 
349       rv = kCHARGEDHADRON ; 
350     break; 
351
352   case -1:    // it's a charged
353     ran = fRan.Rndm() ; 
354     if ( ran <= 0.9996 )
355       rv = kCHARGEDHADRON ; 
356     else 
357       rv = kGAMMA ; 
358
359     break ; 
360   }
361     
362   
363   return rv ;
364 }
365
366 //___________________________________________________________________________
367 void AliPHOSvFast::ResetPoints()
368 {
369   ResetFastRecParticles() ; 
370 }
371
372 //___________________________________________________________________________
373 void AliPHOSvFast::ResetFastRecParticles()
374 {
375   if (fFastRecParticles) 
376     fFastRecParticles->Clear() ;
377   fNRecParticles = 0 ; 
378 }
379
380 //___________________________________________________________________________
381 void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
382 {
383
384   switch (index) {
385   case 0:
386     fBigBoxX = value ; 
387     break ; 
388   case 1:
389     fBigBoxY = value ; 
390     break ; 
391   case 2:
392     fBigBoxZ = value ; 
393     break ; 
394  }
395
396 }
397
398 //____________________________________________________________________________
399 Double_t AliPHOSvFast::SigmaE(Double_t energy)
400 {
401   Double_t rv = -1 ; 
402   
403   rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2) 
404                + TMath::Power(fResPara2/TMath::Sqrt(energy), 2) 
405                + TMath::Power(fResPara3, 2) ) ;  
406
407   return rv * energy ; 
408 }
409
410 //____________________________________________________________________________
411 Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence)
412 {
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
424   Int_t primary =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
425   TLorentzVector lv ; 
426   gMC->TrackPosition(lv) ;
427   TVector3 pos = lv.Vect() ; 
428   Int_t modid  ; 
429   gMC->CurrentVolID(modid);
430   
431   // Makes a reconstructed particle from the primary particle
432
433   TClonesArray * particlelist = gAlice->Particles() ;
434   TParticle * part = (TParticle *)particlelist->At(primary) ;  
435
436   AliPHOSFastRecParticle rp(*part) ;
437
438   // Adds the response of PHOS to the particle
439
440   MakeRecParticle(modid, pos, rp) ;
441   
442   // add the primary particle to the FastRecParticles list
443   
444   AddRecParticle(rp) ;
445
446   // stop the track as soon PHOS is reached
447   
448   gMC->StopTrack() ; 
449
450 }
451