[u/mrichter/AliRoot.git] / PHOS / AliPHOSvFast.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//_________________________________________________________________________
// 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"

// --- Standard library ---

#include <stdio.h>

// --- AliRoot header files ---

#include "AliPHOSvFast.h"
#include "AliRun.h"
#include "AliConst.h"

ClassImp(AliPHOSvFast)

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast()
{
  // ctor

  fFastRecParticles = 0 ; 
  fNRecParticles = 0 ; 
}

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
  AliPHOS(name,title)
{
  // ctor

  // 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->GetCPVBoxSize(1) ) ; 
  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
 
  fFastRecParticles->Delete() ; 
  delete fFastRecParticles ;
  fFastRecParticles = 0 ; 

}

//____________________________________________________________________________
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 ;
  
  new TBRIK( "BigBox", "PHOS box", "void", GetBigBox(0)/2, 
	     GetBigBox(1)/2, 
	     GetBigBox(2)/2 );
  
  // position PHOS into ALICE

  Float_t r = fGeom->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
  Int_t number = 988 ; 
  Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * fGeom->GetIPtoCrystalSurface() ) ) ;
  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 ; 
}

//____________________________________________________________________________
void AliPHOSvFast::CreateGeometry()
{
  // Create the geometry for GEANT
  
  AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
  
  if ( phostmp == NULL ) {
    
    fprintf(stderr, "PHOS detector not found!\n") ;
    return ;
    
  }

  // Get pointer to the array containing media indeces
  Int_t *idtmed = fIdtmed->GetArray() - 699 ;
  
  Float_t bigbox[3] ; 
  bigbox[0] =   GetBigBox(0) / 2.0 ;
  bigbox[1] =   GetBigBox(1) / 2.0 ;
  bigbox[2] =   GetBigBox(2) / 2.0 ;
  
  gMC->Gsvolu("PHOS", "BOX ", idtmed[798], bigbox, 3) ;
  
  // --- Position  PHOS mdules in ALICE setup ---
  
  Int_t idrotm[99] ;
  Double_t const kRADDEG = 180.0 / kPI ;
  
  for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
    
    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->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

}


//____________________________________________________________________________
void AliPHOSvFast::Init(void)
{
  // Prints out an information message
  
  Int_t i;

  printf("\n");
  for(i=0;i<35;i++) printf("*");
  printf(" FAST PHOS_INIT ");
  for(i=0;i<35;i++) printf("*");
  printf("\n");

  // Here the PHOS initialisation code (if any!)

  for(i=0;i<80;i++) printf("*");
  printf("\n");
  
}

//___________________________________________________________________________
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) {
  case 0:
    rv = fBigBoxX ; 
    break ; 
  case 1:
     rv = fBigBoxY ; 
    break ; 
  case 2:
     rv = fBigBoxZ ; 
    break ; 
 }
  return rv ; 
}

//___________________________________________________________________________
void AliPHOSvFast::MakeBranch(Option_t* opt, char *file)
{  
  // Create new branch in the current reconstructed Root Tree
 
  AliDetector::MakeBranch(opt,file) ;
  
  char branchname[10];
  sprintf(branchname,"%s",GetName());
  char *cd = strstr(opt,"R");
  
  if (fFastRecParticles && gAlice->TreeR() && cd) {
    gAlice->MakeBranchInTree(gAlice->TreeR(), 
                             branchname, &fFastRecParticles, fBufferSize, file);
  }
}

//____________________________________________________________________________
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::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 ; 
    break ; 
  case 1:
    fBigBoxY = value ; 
    break ; 
  case 2:
    fBigBoxZ = value ; 
    break ; 
 }

}

//____________________________________________________________________________
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() ); 
  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) ;  

  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
  
  gMC->StopTrack() ; 

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