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