aa8d558902d4089bacf7755b3d909d728f6a21ef
[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   // The branch of fHits is automatically split
224   //
225   AliDetector::MakeBranch(opt) ;
226   
227   char branchname[10];
228   sprintf(branchname,"%s",GetName());
229   char *cd = strstr(opt,"R");
230   
231   if (fFastRecParticles && gAlice->TreeR() && cd) {
232     gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
233   }
234 }
235
236 //____________________________________________________________________________
237 Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
238 {  
239   Double_t sigma  = SigmaE(energy) ; 
240   return  fRan.Gaus(energy, sigma) ;   
241 }
242
243 //____________________________________________________________________________
244 TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi)
245 {
246   TVector3 newpos ;
247   Double_t sigma = SigmaP( energy, theta*180./TMath::Pi() ) ;
248   Double_t x = fRan.Gaus( pos.X(), sigma ) ;
249   sigma = SigmaP( energy, phi*180./TMath::Pi() ) ;
250   Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
251   Double_t y = pos.Y() ; 
252   
253   newpos.SetX(x) ; 
254   newpos.SetY(y) ; 
255   newpos.SetZ(z) ; 
256               
257   return newpos ; 
258 }
259
260 //____________________________________________________________________________
261 void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
262 {
263   // get the detected type of particle
264  
265   Int_t type = MakeType( rp ) ;
266   rp.SetType(type) ;
267
268   
269   // get the detected energy
270
271   TLorentzVector momentum ;  
272   rp.Momentum(momentum) ; 
273   Double_t kineticenergy = TMath::Sqrt( TMath::Power(momentum.E(), 2) - TMath::Power(rp.GetMass(), 2) ) ; 
274   Double_t modifiedkineticenergy = MakeEnergy(kineticenergy ) ;
275   Double_t modifiedenergy = TMath::Sqrt( TMath::Power(modifiedkineticenergy, 2)  
276                                          + TMath::Power( rp.GetMass(), 2) ) ;
277  
278   // get the angle of incidence 
279   
280   Double_t incidencetheta = 90. * TMath::Pi() /180 - rp.Theta() ; 
281   Double_t incidencephi   = ( 270 + fGeom->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;   
282
283   // get the detected direction
284   
285   TVector3 modifiedposition = MakePosition(kineticenergy, pos, incidencetheta, incidencephi) ; 
286   modifiedposition *= modifiedkineticenergy / modifiedposition.Mag() ; 
287
288   // Set the modified 4-momentum of the reconstructed particle
289
290   rp.SetMomentum(modifiedposition.X(), modifiedposition.Y(), modifiedposition.Z(), modifiedenergy) ; 
291
292  }
293
294 //____________________________________________________________________________
295 Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
296 {
297   Int_t rv =  kUNDEFINED ;
298   Int_t charge = (Int_t)rp.GetPDG()->Charge() ;
299   Int_t test ; 
300   Float_t ran ; 
301   if ( charge == 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) ) 
302     test = - 1 ;
303   else
304     test = rp.GetPdgCode() ; 
305
306   switch (test) { 
307
308   case 22:    // it's a photon
309     ran = fRan.Rndm() ; 
310     if ( ran <= 0.5 )  // 50 % 
311       rv = kGAMMA ; 
312     else {
313       ran = fRan.Rndm() ; 
314       if( ran <= 0.9498 )
315         rv = kNEUTRALEM ; 
316       else
317         rv = kNEUTRALHADRON ; 
318     }     
319     break ; 
320
321   case 2112:  // it's a neutron
322     ran = fRan.Rndm() ; 
323     if ( ran <= 0.9998 )
324       rv = kNEUTRALHADRON ; 
325     else 
326       rv = kNEUTRALEM ; 
327     break ; 
328
329   case -2112: // it's a anti-neutron
330     ran = fRan.Rndm() ; 
331     if ( ran <= 0.9984 )
332       rv = kNEUTRALHADRON ; 
333     else 
334       rv = kNEUTRALEM ; 
335     break ; 
336     
337   case 11:    // it's a electron
338     ran = fRan.Rndm() ; 
339     if ( ran <= 0.9996 )
340       rv = kELECTRON ; 
341     else 
342       rv = kCHARGEDHADRON ; 
343     break; 
344
345   case -11:   // it's a electron
346     ran = fRan.Rndm() ; 
347     if ( ran <= 0.9996 )
348       rv = kELECTRON ; 
349     else 
350       rv = kCHARGEDHADRON ; 
351     break; 
352
353   case -1:    // it's a charged
354     ran = fRan.Rndm() ; 
355     if ( ran <= 0.9996 )
356       rv = kCHARGEDHADRON ; 
357     else 
358       rv = kGAMMA ; 
359
360     break ; 
361   }
362     
363   
364   return rv ;
365 }
366
367 //___________________________________________________________________________
368 void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
369 {
370
371   switch (index) {
372   case 0:
373     fBigBoxX = value ; 
374     break ; 
375   case 1:
376     fBigBoxY = value ; 
377     break ; 
378   case 2:
379     fBigBoxZ = value ; 
380     break ; 
381  }
382
383 }
384
385 //____________________________________________________________________________
386 Double_t AliPHOSvFast::SigmaE(Double_t energy)
387 {
388   Double_t rv = -1 ; 
389   
390   rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2) 
391                + TMath::Power(fResPara2/TMath::Sqrt(energy), 2) 
392                + TMath::Power(fResPara3, 2) ) ;  
393
394   return rv * energy ; 
395 }
396
397 //____________________________________________________________________________
398 Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence)
399 {
400  
401   Double_t paraA = fPosParaA0 + fPosParaA1 * incidence ; 
402   Double_t paraB = fPosParaB0 + fPosParaB1 * incidence + fPosParaB2 * incidence * incidence ; 
403
404   return ( paraA / TMath::Sqrt(energy) + paraB ) * 0.1   ; // in cm  
405 }
406
407 //____________________________________________________________________________
408 void AliPHOSvFast::StepManager(void)
409 {
410
411   Int_t primary =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
412   TLorentzVector lv ; 
413   gMC->TrackPosition(lv) ;
414   TVector3 pos = lv.Vect() ; 
415   Int_t modid  ; 
416   gMC->CurrentVolID(modid);
417   
418   // Makes a reconstructed particle from the primary particle
419
420   TClonesArray * particlelist = gAlice->Particles() ;
421   TParticle * part = (TParticle *)particlelist->At(primary) ;  
422
423   AliPHOSFastRecParticle rp(*part) ;
424
425   // Adds the response of PHOS to the particle
426
427   MakeRecParticle(modid, pos, rp) ;
428   
429   // add the primary particle to the FastRecParticles list
430   
431   AddRecParticle(rp) ;
432
433   // stop the track as soon PHOS is reached
434   
435   gMC->StopTrack() ; 
436
437 }
438