Put back the two bellows in front of the absorber.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSvFast.cxx
index 08c5216..8add01e 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //_________________________________________________________________________
-// 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 ---
 
 
 // --- Standard library ---
 
-#include <cstdio>
-#include <cassert>
+#include <stdio.h>
 
 // --- AliRoot header files ---
 
 #include "AliPHOSvFast.h"
-#include "AliPHOSReconstructioner.h"
 #include "AliRun.h"
 #include "AliConst.h"
 
@@ -41,7 +45,9 @@ ClassImp(AliPHOSvFast)
 //____________________________________________________________________________
 AliPHOSvFast::AliPHOSvFast()
 {
-  fRecParticles = 0 ; 
+  // ctor
+
+  fFastRecParticles = 0 ; 
   fNRecParticles = 0 ; 
 }
 
@@ -49,6 +55,8 @@ AliPHOSvFast::AliPHOSvFast()
 AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
   AliPHOS(name,title)
 {
+  // ctor
+
   // gets an instance of the geometry parameters class  
    
   fGeom =  AliPHOSGeometry::GetInstance(title, "") ; 
@@ -63,32 +71,61 @@ AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
   SetBigBox(2, fGeom->GetOuterBoxSize(0) ); 
 
   fNRecParticles = 0 ; 
+  fFastRecParticles = new AliPHOSFastRecParticle::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 ; 
 }
 
 //____________________________________________________________________________
 AliPHOSvFast::~AliPHOSvFast()
 {
+  // dtor
  
-  fRecParticles->Delete() ; 
-  delete fRecParticles ;
-  fRecParticles = 0 ; 
+  fFastRecParticles->Delete() ; 
+  delete fFastRecParticles ;
+  fFastRecParticles = 0 ; 
 
 }
 
 //____________________________________________________________________________
-void AliPHOSvFast::AddRecParticle(Int_t primary)
-{
-   TClonesArray * particlelist = gAlice->Particles() ;
-   TParticle * part = (TParticle *)particlelist->At(primary) ;  
-   cout <<  " AliPHOSvFast::AddRecParticle " << part->GetName() << endl ; 
+void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
+{  
+  // Add a virtually reconstructed particle to the list 
+
+  new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
+  fNRecParticles++ ; 
 }
 
 //____________________________________________________________________________
 void AliPHOSvFast::BuildGeometry()
 {
-
   // Build the PHOS geometry for the ROOT display
-  
+   //BEGIN_HTML
+  /*
+    <H2>
+     PHOS FAST in ALICE displayed by root
+    </H2>
+    <H4> All Views </H4>
+    <P>
+    <CENTER>
+    <IMG Align=BOTTOM ALT="Fast All Views" SRC="../images/AliPHOSvFastAllViews.gif"> 
+    </CENTER></P>
+    <H4> Front View </H4>
+    <P>
+    <CENTER>
+    <IMG Align=BOTTOM ALT="Fast Front View" SRC="../images/AliPHOSvFastFrontView.gif"> 
+    </CENTER></P>
+  */
+  //END_HTML  
   const Int_t kColorPHOS = kRed ;
   
   Double_t const kRADDEG = 180.0 / kPI ;
@@ -99,9 +136,9 @@ void AliPHOSvFast::BuildGeometry()
   
   // position PHOS into ALICE
 
-  Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
+  Float_t r = fGeom->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
   Int_t number = 988 ; 
-  Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * fGeom->GetIPtoOuterCoverDistance() ) ) ;
+  Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * fGeom->GetIPtoCrystalSurface() ) ) ;
   pphi *= kRADDEG ;
   TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
  
@@ -127,9 +164,10 @@ void AliPHOSvFast::BuildGeometry()
 //____________________________________________________________________________
 void AliPHOSvFast::CreateGeometry()
 {
-
+  // Create the geometry for GEANT
+  
   AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
-
+  
   if ( phostmp == NULL ) {
     
     fprintf(stderr, "PHOS detector not found!\n") ;
@@ -157,11 +195,10 @@ void AliPHOSvFast::CreateGeometry()
     Float_t angle = fGeom->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 = fGeom->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
 
     Float_t xP1 = r * TMath::Sin( angle / kRADDEG ) ;
     Float_t yP1 = -r * TMath::Cos( angle / kRADDEG ) ;
-
     gMC->Gspos("PHOS", i, "ALIC", xP1, yP1, 0.0, idrotm[i-1], "ONLY") ;
  
   } // for GetNModules
@@ -172,7 +209,8 @@ void AliPHOSvFast::CreateGeometry()
 //____________________________________________________________________________
 void AliPHOSvFast::Init(void)
 {
+  // Prints out an information message
+  
   Int_t i;
 
   printf("\n");
@@ -191,6 +229,8 @@ void AliPHOSvFast::Init(void)
 //___________________________________________________________________________
 Float_t AliPHOSvFast::GetBigBox(Int_t index)
 {
+  // Get the X, Y or Z dimension of the box describing a PHOS module
+  
   Float_t rv = 0 ; 
 
   switch (index) {
@@ -210,26 +250,183 @@ Float_t AliPHOSvFast::GetBigBox(Int_t index)
 //___________________________________________________________________________
 void AliPHOSvFast::MakeBranch(Option_t* opt)
 {  
-  //
-  // Create a new branch in the current Root Tree
-  // The branch of fHits is automatically split
-  //
+  // Create new branch in the current reconstructed Root Tree
   AliDetector::MakeBranch(opt) ;
   
   char branchname[10];
   sprintf(branchname,"%s",GetName());
-  char *cd = strstr(opt,"D");
+  char *cd = strstr(opt,"R");
   
-  if (fDigits && gAlice->TreeD() && cd) {
-    gAlice->TreeD()->Branch(branchname, &fRecParticles, fBufferSize);
-    //    printf("* AliPHOS::MakeBranch * Making Branch %s for RecParticles \n",branchname);
+  if (fFastRecParticles && gAlice->TreeR() && cd) {
+    gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
   }
 }
 
+//____________________________________________________________________________
+Double_t AliPHOSvFast::MakeEnergy(const 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)
+{
+  // 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 ) ;
+  sigma = SigmaP( energy, phi*180./TMath::Pi() ) ;
+  Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
+  Double_t y = pos.Y() ; 
+  
+  newpos.SetX(x) ; 
+  newpos.SetY(y) ; 
+  newpos.SetZ(z) ; 
+             
+  return newpos ; 
+}
+
+//____________________________________________________________________________
+void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
+{
+  // 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) ;
+
+  
+  // get the detected energy
+
+  TLorentzVector momentum ;  
+  rp.Momentum(momentum) ; 
+  Double_t kineticenergy = TMath::Sqrt( TMath::Power(momentum.E(), 2) - TMath::Power(rp.GetMass(), 2) ) ; 
+  Double_t modifiedkineticenergy = MakeEnergy(kineticenergy ) ;
+  Double_t modifiedenergy = TMath::Sqrt( TMath::Power(modifiedkineticenergy, 2)  
+                                        + TMath::Power( rp.GetMass(), 2) ) ;
+  // 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() ;   
+
+  // get the detected direction
+  
+  TVector3 modifiedposition = MakePosition(kineticenergy, pos, incidencetheta, incidencephi) ; 
+  modifiedposition *= modifiedkineticenergy / modifiedposition.Mag() ; 
+
+  // Set the modified 4-momentum of the reconstructed particle
+
+  rp.SetMomentum(modifiedposition.X(), modifiedposition.Y(), modifiedposition.Z(), modifiedenergy) ; 
+
+ }
+
+//____________________________________________________________________________
+Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
+{
+  // 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 ) ) 
+    test = - 1 ;
+  else
+    test = rp.GetPdgCode() ; 
+
+  switch (test) { 
+
+  case 22:    // it's a photon
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.5 )  // 50 % 
+      rv =  AliPHOSFastRecParticle::kGAMMA ; 
+    else {
+      ran = fRan.Rndm() ; 
+      if( ran <= 0.9498 )
+       rv =  AliPHOSFastRecParticle::kNEUTRALEM ; 
+      else
+       rv =  AliPHOSFastRecParticle::kNEUTRALHA ; 
+    }     
+    break ; 
+
+  case 2112:  // it's a neutron
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.9998 )
+      rv =  AliPHOSFastRecParticle::kNEUTRALHA ; 
+    else 
+      rv = AliPHOSFastRecParticle::kNEUTRALEM ; 
+    break ; 
+
+  case -2112: // it's a anti-neutron
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.9984 )
+      rv =  AliPHOSFastRecParticle::kNEUTRALHA ; 
+    else 
+      rv =  AliPHOSFastRecParticle::kNEUTRALEM ; 
+    break ; 
+    
+  case 11:    // it's a electron
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.9996 )
+      rv =  AliPHOSFastRecParticle::kELECTRON ; 
+    else 
+      rv =  AliPHOSFastRecParticle::kCHARGEDHA ; 
+    break; 
+
+  case -11:   // it's a positon
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.9996 )
+      rv =  AliPHOSFastRecParticle::kELECTRON ; 
+    else 
+      rv =  AliPHOSFastRecParticle::kCHARGEDHA ; 
+    break; 
+
+  case -1:    // it's a charged
+    ran = fRan.Rndm() ; 
+    if ( ran <= 0.9996 )
+      rv =  AliPHOSFastRecParticle::kCHARGEDHA ; 
+    else 
+      rv =  AliPHOSFastRecParticle::kGAMMA ; 
+
+    break ; 
+  }
+    
+  
+  return rv ;
+}
+
 //___________________________________________________________________________
-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 ; 
@@ -245,16 +442,60 @@ void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
 }
 
 //____________________________________________________________________________
-void AliPHOSvFast::StepManager(void)
+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) 
+              + TMath::Power(fResPara2/TMath::Sqrt(energy), 2) 
+              + TMath::Power(fResPara3, 2) ) ;  
+
+  return rv * energy ; 
+}
+
+//____________________________________________________________________________
+Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence)
 {
+  // Calculates the energy dependent position resolution 
+
+  Double_t paraA = fPosParaA0 + fPosParaA1 * incidence ; 
+  Double_t paraB = fPosParaB0 + fPosParaB1 * incidence + fPosParaB2 * incidence * incidence ; 
+
+  return ( paraA / TMath::Sqrt(energy) + paraB ) * 0.1   ; // in cm  
+}
 
+//____________________________________________________________________________
+void AliPHOSvFast::StepManager(void)
+{
+  // Only verifies if the particle reaches PHOS and stops the tracking 
+  
   Int_t primary =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
-  TString name = fGeom->GetName() ; 
+  TLorentzVector lv ; 
+  gMC->TrackPosition(lv) ;
+  TVector3 pos = lv.Vect() ; 
+  Int_t modid  ; 
+  gMC->CurrentVolID(modid);
+  
+  // Makes a reconstructed particle from the primary particle
+
+  TClonesArray * particlelist = gAlice->Particles() ;
+  TParticle * part = (TParticle *)particlelist->At(primary) ;  
 
-  // add the primary particle to the RecParticles list
+  AliPHOSFastRecParticle rp(*part) ;
+  rp.SetPrimary(primary) ; 
+
+  // Adds the response of PHOS to the particle
+
+  MakeRecParticle(modid, pos, rp) ;
+  
+  // add the primary particle to the FastRecParticles list
+  
+  AddRecParticle(rp) ;
+
+  // stop the track as soon PHOS is reached
   
-  AddRecParticle(primary);
-  fNRecParticles++ ; 
   gMC->StopTrack() ; 
 
 }