Updates for the ROOT trunk
[u/mrichter/AliRoot.git] / PHOS / AliPHOSvFast.cxx
index aa8d558..32ae66f 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+  /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.30  2006/09/13 07:31:01  kharlov
+ * Effective C++ corrections (T.Pocheptsov)
+ *
+ * Revision 1.29  2005/05/28 14:19:05  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
-// Manager class for PHOS version for fast simulations
-//*-- Author : Y. Schutz SUBATECH 
-//////////////////////////////////////////////////////////////////////////////
+// Implementation of the PHOS manager class for fast simulations     
+// Tracks particles until the reach a grossly designed PHOS module
+// Modify the particles property (momentum, energy, type) according to
+//  the PHOS response function. The result is called a virtual reconstructed
+//  particle.                
+//
+//*-- Author: Yves Schutz (SUBATECH)
 
 // --- ROOT system ---
-
-#include "TBRIK.h"
-#include "TNode.h"
-#include "TParticle.h"
+#include <TGeometry.h>
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include "TClonesArray.h" 
+#include <TVirtualMC.h>
 
 // --- Standard library ---
 
-#include <cstdio>
-#include <cassert>
-
 // --- AliRoot header files ---
-
+#include "AliPHOSFastRecParticle.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSLoader.h"
 #include "AliPHOSvFast.h"
 #include "AliRun.h"
-#include "AliConst.h"
 
 ClassImp(AliPHOSvFast)
 
-//____________________________________________________________________________
-AliPHOSvFast::AliPHOSvFast()
+AliPHOSvFast::AliPHOSvFast() :
+  fBigBoxX(0.),
+  fBigBoxY(0.),
+  fBigBoxZ(0.),
+  fFastRecParticles(0),
+  fNRecParticles(0),
+  fRan(0),
+  fResPara1(0.),
+  fResPara2(0.),
+  fResPara3(0.),
+  fPosParaA0(0.),
+  fPosParaA1(0.),
+  fPosParaB0(0.),
+  fPosParaB1(0.),
+  fPosParaB2(0.)    
 {
-  fFastRecParticles = 0 ; 
-  fNRecParticles = 0 ; 
+  // default ctor : initialize data member
 }
 
 //____________________________________________________________________________
 AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
-  AliPHOS(name,title)
+  AliPHOS(name,title),
+  fBigBoxX(0.),
+  fBigBoxY(0.),
+  fBigBoxZ(0.),
+  fFastRecParticles(new AliPHOSFastRecParticle::FastRecParticlesList("AliPHOSFastRecParticle", 100)),
+  fNRecParticles(0),
+  fRan(0),
+  fResPara1(0.030), // GeV
+  fResPara2(0.00003),
+  fResPara3(0.00001),
+  fPosParaA0(2.87),    // mm
+  fPosParaA1(-0.0975),
+  fPosParaB0(0.257),
+  fPosParaB1(0.137),
+  fPosParaB2(0.00619)
 {
-  // gets an instance of the geometry parameters class  
-   
-  fGeom =  AliPHOSGeometry::GetInstance(title, "") ; 
-
-  if (fGeom->IsInitialized() ) 
-    cout << "AliPHOSvFast : PHOS geometry intialized for " << fGeom->GetName() << endl ;
-  else
-    cout << "AliPHOSvFast : PHOS geometry initialization failed !" << endl ;   
-  
-  SetBigBox(0, fGeom->GetOuterBoxSize(0) ) ;
-  SetBigBox(1, fGeom->GetOuterBoxSize(1) + fGeom->GetPPSDBoxSize(1) ) ; 
-  SetBigBox(2, fGeom->GetOuterBoxSize(0) ); 
-
-  fNRecParticles = 0 ; 
-  fFastRecParticles = new FastRecParticlesList("AliPHOSFastRecParticle", 100) ;
-
-  fResPara1 = 0.030 ;    // GeV
-  fResPara2 = 0.00003 ; 
-  fResPara3 = 0.00001 ; 
-
-  fPosParaA0 = 2.87 ;    // mm
-  fPosParaA1 = -0.0975 ;  
-  fPosParaB0 = 0.257 ;   
-  fPosParaB1 = 0.137 ; 
-  fPosParaB2 = 0.00619 ; 
+  // ctor
+  // create the Loader 
+  SetBigBox(0, GetGeometry()->GetOuterBoxSize(0) ) ;
+  SetBigBox(1, GetGeometry()->GetOuterBoxSize(3) + GetGeometry()->GetCPVBoxSize(1) ) ; 
+  SetBigBox(2, GetGeometry()->GetOuterBoxSize(2) ); 
 }
 
 //____________________________________________________________________________
 AliPHOSvFast::~AliPHOSvFast()
 {
+  // dtor
  
   fFastRecParticles->Delete() ; 
   delete fFastRecParticles ;
@@ -88,57 +112,19 @@ AliPHOSvFast::~AliPHOSvFast()
 //____________________________________________________________________________
 void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
 {  
-   new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
-   fNRecParticles++ ; 
-}
+  // Add a virtually reconstructed particle to the list 
 
-//____________________________________________________________________________
-void AliPHOSvFast::BuildGeometry()
-{
-
-  // Build the PHOS geometry for the ROOT display
-  
-  const Int_t kColorPHOS = kRed ;
-  
-  Double_t const kRADDEG = 180.0 / kPI ;
-  
-  new TBRIK( "BigBox", "PHOS box", "void", GetBigBox(0)/2, 
-            GetBigBox(1)/2, 
-            GetBigBox(2)/2 );
-  
-  // position PHOS into ALICE
-
-  Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
-  Int_t number = 988 ; 
-  Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * fGeom->GetIPtoOuterCoverDistance() ) ) ;
-  pphi *= kRADDEG ;
-  TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
-  char * nodename = new char[20] ;  
-  char * rotname  = new char[20] ; 
-
-  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) { 
-   Float_t angle = pphi * 2 * ( i - fGeom->GetNModules() / 2.0 - 0.5 ) ;
-   sprintf(rotname, "%s%d", "rot", number++) ;
-   new TRotMatrix(rotname, rotname, 90, angle, 90, 90 + angle, 0, 0);
-   top->cd();
-   sprintf(nodename,"%s%d", "Module", i) ;    
-   Float_t x =  r * TMath::Sin( angle / kRADDEG ) ;
-   Float_t y = -r * TMath::Cos( angle / kRADDEG ) ;
-   TNode * bigboxnode = new TNode(nodename, nodename, "BigBox", x, y, 0, rotname ) ;
-   bigboxnode->SetLineColor(kColorPHOS) ;
-   fNodes->Add(bigboxnode) ;
-  }
-  delete[] nodename ; 
-  delete[] rotname ; 
+  new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
+  fNRecParticles++ ; 
 }
 
 //____________________________________________________________________________
 void AliPHOSvFast::CreateGeometry()
 {
-
+  // Create the geometry for GEANT
+  
   AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
-
+  
   if ( phostmp == NULL ) {
     
     fprintf(stderr, "PHOS detector not found!\n") ;
@@ -159,14 +145,14 @@ void AliPHOSvFast::CreateGeometry()
   // --- Position  PHOS mdules in ALICE setup ---
   
   Int_t idrotm[99] ;
-  Double_t const kRADDEG = 180.0 / kPI ;
+  Double_t const kRADDEG = 180.0 / TMath::Pi() ;
   
-  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
+  for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
     
-    Float_t angle = fGeom->GetPHOSAngle(i) ;
+    Float_t angle = GetGeometry()->GetPHOSAngle(i) ;
     AliMatrix(idrotm[i-1], 90.0, angle, 90.0, 90.0+angle, 0.0, 0.0) ;
  
-    Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
+    Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
 
     Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
     Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
@@ -180,7 +166,8 @@ void AliPHOSvFast::CreateGeometry()
 //____________________________________________________________________________
 void AliPHOSvFast::Init(void)
 {
+  // Prints out an information message
+  
   Int_t i;
 
   printf("\n");
@@ -197,8 +184,10 @@ void AliPHOSvFast::Init(void)
 }
 
 //___________________________________________________________________________
-Float_t AliPHOSvFast::GetBigBox(Int_t index)
+Float_t AliPHOSvFast::GetBigBox(Int_t index) const
 {
+  // Get the X, Y or Z dimension of the box describing a PHOS module
+  
   Float_t rv = 0 ; 
 
   switch (index) {
@@ -214,35 +203,35 @@ Float_t AliPHOSvFast::GetBigBox(Int_t index)
  }
   return rv ; 
 }
-
 //___________________________________________________________________________
+
 void AliPHOSvFast::MakeBranch(Option_t* opt)
 {  
-  //
-  // Create a new branch in the current Root Tree
-  // The branch of fHits is automatically split
-  //
-  AliDetector::MakeBranch(opt) ;
-  
-  char branchname[10];
-  sprintf(branchname,"%s",GetName());
-  char *cd = strstr(opt,"R");
+  // Create new branch in the current reconstructed Root Tree
+  AliDetector::MakeBranch(opt);
+  const char *cd = strstr(opt,"R");
   
-  if (fFastRecParticles && gAlice->TreeR() && cd) {
-    gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
+  if (fFastRecParticles && fLoader->TreeR() && cd) {
+    MakeBranchInTree(fLoader->TreeR(), GetName(), &fFastRecParticles, fBufferSize, 0);
   }
 }
-
 //____________________________________________________________________________
-Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
+
+Double_t AliPHOSvFast::MakeEnergy(Double_t energy)
 {  
+  // Smears the energy according to the energy dependent energy resolution.
+  // A gaussian distribution is assumed
+
   Double_t sigma  = SigmaE(energy) ; 
   return  fRan.Gaus(energy, sigma) ;   
 }
-
 //____________________________________________________________________________
-TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi)
+
+TVector3 AliPHOSvFast::MakePosition(Double_t energy, TVector3 pos, Double_t theta, Double_t phi)
 {
+  // Smears the impact position according to the energy dependent position resolution
+  // A gaussian position distribution is assumed
+
   TVector3 newpos ;
   Double_t sigma = SigmaP( energy, theta*180./TMath::Pi() ) ;
   Double_t x = fRan.Gaus( pos.X(), sigma ) ;
@@ -258,10 +247,12 @@ TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, c
 }
 
 //____________________________________________________________________________
-void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
+void AliPHOSvFast::MakeRecParticle(Int_t modid, TVector3 pos, AliPHOSFastRecParticle & rp)
 {
-  // get the detected type of particle
+  // Modify the primary particle properties according
+  //  1. the response function of PHOS
+  //  2. the performance of the EMC+PPSD setup
+  
   Int_t type = MakeType( rp ) ;
   rp.SetType(type) ;
 
@@ -278,7 +269,7 @@ void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHO
   // get the angle of incidence 
   
   Double_t incidencetheta = 90. * TMath::Pi() /180 - rp.Theta() ; 
-  Double_t incidencephi   = ( 270 + fGeom->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;   
+  Double_t incidencephi   = ( 270 + GetGeometry()->GetPHOSAngle(modid) ) * TMath::Pi() / 180. - rp.Phi() ;   
 
   // get the detected direction
   
@@ -294,68 +285,67 @@ void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHO
 //____________________________________________________________________________
 Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
 {
-  Int_t rv =  kUNDEFINED ;
+  // Generate a particle type using the performance of the EMC+PPSD setup
+
+  Int_t rv =   AliPHOSFastRecParticle::kUNDEFINED ;
   Int_t charge = (Int_t)rp.GetPDG()->Charge() ;
   Int_t test ; 
   Float_t ran ; 
-  if ( charge == 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) ) 
+  if ( charge != 0 && ( TMath::Abs(rp.GetPdgCode()) != 11 ) ) 
     test = - 1 ;
   else
     test = rp.GetPdgCode() ; 
 
+  Fatal("MakeType", "SHOULD NOT BE USED until values of probabilities are properly set ") ;
+   // NB: ALL VALUES SHOULD BE CHECKED !!!!
   switch (test) { 
 
-  case 22:    // it's a photon
+  case 22:    // it's a photon              // NB: ALL VALUES SHOLD BE CHECKED !!!!
     ran = fRan.Rndm() ; 
-    if ( ran <= 0.5 )  // 50 % 
-      rv = kGAMMA ; 
-    else {
-      ran = fRan.Rndm() ; 
-      if( ran <= 0.9498 )
-       rv = kNEUTRALEM ; 
-      else
-       rv = kNEUTRALHADRON ; 
-    }     
+    if( ran <= 0.9498 )
+      rv =  AliPHOSFastRecParticle::kNEUTRALHAFAST ; 
+    else
+      rv =  AliPHOSFastRecParticle::kNEUTRALEMFAST ;     
     break ; 
 
   case 2112:  // it's a neutron
     ran = fRan.Rndm() ; 
     if ( ran <= 0.9998 )
-      rv = kNEUTRALHADRON ; 
+      rv =  AliPHOSFastRecParticle::kNEUTRALHASLOW ; 
     else 
-      rv = kNEUTRALEM ; 
+      rv = AliPHOSFastRecParticle::kNEUTRALEMSLOW ; 
     break ; 
-
+    
   case -2112: // it's a anti-neutron
     ran = fRan.Rndm() ; 
     if ( ran <= 0.9984 )
-      rv = kNEUTRALHADRON ; 
+      rv =  AliPHOSFastRecParticle::kNEUTRALHASLOW ; 
     else 
-      rv = kNEUTRALEM ; 
+      rv =  AliPHOSFastRecParticle::kNEUTRALEMSLOW ; 
     break ; 
     
   case 11:    // it's a electron
     ran = fRan.Rndm() ; 
     if ( ran <= 0.9996 )
-      rv = kELECTRON ; 
+      rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ; 
     else 
-      rv = kCHARGEDHADRON ; 
+      rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ; 
     break; 
 
-  case -11:   // it's a electron
+  case -11:   // it's a positon
     ran = fRan.Rndm() ; 
     if ( ran <= 0.9996 )
-      rv = kELECTRON ; 
+      rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ; 
     else 
-      rv = kCHARGEDHADRON ; 
+      rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ; 
     break; 
 
   case -1:    // it's a charged
     ran = fRan.Rndm() ; 
     if ( ran <= 0.9996 )
-      rv = kCHARGEDHADRON ; 
+      rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ; 
     else 
-      rv = kGAMMA ; 
+      rv =  AliPHOSFastRecParticle::kNEUTRALHAFAST ; 
 
     break ; 
   }
@@ -365,9 +355,28 @@ Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
 }
 
 //___________________________________________________________________________
-void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
+void AliPHOSvFast::ResetPoints()
 {
+  // This overloads the method in AliDetector
+  
+  ResetFastRecParticles() ; 
+}
 
+//___________________________________________________________________________
+void AliPHOSvFast::ResetFastRecParticles()
+{
+  // Resets the list of virtual reconstructed particles
+  if (fFastRecParticles) 
+    fFastRecParticles->Clear() ;
+  fNRecParticles = 0 ; 
+}
+
+//___________________________________________________________________________
+void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
+{
+  // Set the size of the Box describing a PHOS module
+  
   switch (index) {
   case 0:
     fBigBoxX = value ; 
@@ -385,6 +394,8 @@ void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
 //____________________________________________________________________________
 Double_t AliPHOSvFast::SigmaE(Double_t energy)
 {
+  // Calculates the energy dependent energy resolution
+  
   Double_t rv = -1 ; 
   
   rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2) 
@@ -395,9 +406,10 @@ Double_t AliPHOSvFast::SigmaE(Double_t energy)
 }
 
 //____________________________________________________________________________
-Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence)
+Double_t AliPHOSvFast::SigmaP(Double_t energy, Double_t incidence)
 {
+  // Calculates the energy dependent position resolution 
+
   Double_t paraA = fPosParaA0 + fPosParaA1 * incidence ; 
   Double_t paraB = fPosParaB0 + fPosParaB1 * incidence + fPosParaB2 * incidence * incidence ; 
 
@@ -407,29 +419,38 @@ Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence)
 //____________________________________________________________________________
 void AliPHOSvFast::StepManager(void)
 {
+  // Only verifies if the particle reaches PHOS and stops the tracking 
 
-  Int_t primary =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
   TLorentzVector lv ; 
   gMC->TrackPosition(lv) ;
   TVector3 pos = lv.Vect() ; 
   Int_t modid  ; 
   gMC->CurrentVolID(modid);
   
-  // Makes a reconstructed particle from the primary particle
+  Float_t energy = gMC->Etot() ; //Total energy of current track
 
-  TClonesArray * particlelist = gAlice->Particles() ;
-  TParticle * part = (TParticle *)particlelist->At(primary) ;  
+  //Calculating mass of current particle
+  TDatabasePDG * pdg = TDatabasePDG::Instance() ;
+  TParticlePDG * partPDG = pdg->GetParticle(gMC->TrackPid()) ;
+  Float_t mass = partPDG->Mass() ;
 
-  AliPHOSFastRecParticle rp(*part) ;
+  if(energy > mass){
+    pos.SetMag(TMath::Sqrt(energy*energy-mass*mass)) ;
+    TLorentzVector pTrack(pos, energy) ;  
 
-  // Adds the response of PHOS to the particle
+    TParticle * part = new TParticle(gMC->TrackPid(), 0,-1,-1,-1,-1, pTrack, lv)  ;
+        
+    AliPHOSFastRecParticle rp(*part) ;
 
-  MakeRecParticle(modid, pos, rp) ;
-  
-  // add the primary particle to the FastRecParticles list
+    // Adds the response of PHOS to the particle
+    MakeRecParticle(modid, pos, rp) ;
+    
+    // add the `track' particle to the FastRecParticles list
   
-  AddRecParticle(rp) ;
+    AddRecParticle(rp) ;
 
+    part->Delete() ;
+  }
   // stop the track as soon PHOS is reached
   
   gMC->StopTrack() ;