[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.                  *
 **************************************************************************/

//_________________________________________________________________________
// Manager class for PHOS version for fast simulations
//*-- Author : Y. Schutz SUBATECH 
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---

#include "TBRIK.h"
#include "TNode.h"
#include "TParticle.h"

// --- Standard library ---

#include <cstdio>
#include <cassert>

// --- AliRoot header files ---

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

ClassImp(AliPHOSvFast)

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast()
{
  fFastRecParticles = 0 ; 
  fNRecParticles = 0 ; 
}

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
  AliPHOS(name,title)
{
  // 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 ; 
}

//____________________________________________________________________________
AliPHOSvFast::~AliPHOSvFast()
{
 
  fFastRecParticles->Delete() ; 
  delete fFastRecParticles ;
  fFastRecParticles = 0 ; 

}

//____________________________________________________________________________
void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
{  
   new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
   fNRecParticles++ ; 
}

//____________________________________________________________________________
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 ; 
}

//____________________________________________________________________________
void AliPHOSvFast::CreateGeometry()
{

  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->GetIPtoOuterCoverDistance() + 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)
{
 
  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)
{
  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)
{  
  //
  // 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");
  
  if (fFastRecParticles && gAlice->TreeR() && cd) {
    gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
  }
}

//____________________________________________________________________________
Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
{  
  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 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)
{
  // get the detected type of particle
 
  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 )
{
  Int_t rv =  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 = kGAMMA ; 
    else {
      ran = fRan.Rndm() ; 
      if( ran <= 0.9498 )
	rv = kNEUTRALEM ; 
      else
	rv = kNEUTRALHADRON ; 
    }     
    break ; 

  case 2112:  // it's a neutron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9998 )
      rv = kNEUTRALHADRON ; 
    else 
      rv = kNEUTRALEM ; 
    break ; 

  case -2112: // it's a anti-neutron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9984 )
      rv = kNEUTRALHADRON ; 
    else 
      rv = kNEUTRALEM ; 
    break ; 
    
  case 11:    // it's a electron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9996 )
      rv = kELECTRON ; 
    else 
      rv = kCHARGEDHADRON ; 
    break; 

  case -11:   // it's a electron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9996 )
      rv = kELECTRON ; 
    else 
      rv = kCHARGEDHADRON ; 
    break; 

  case -1:    // it's a charged
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9996 )
      rv = kCHARGEDHADRON ; 
    else 
      rv = kGAMMA ; 

    break ; 
  }
    
  
  return rv ;
}

//___________________________________________________________________________
void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
{

  switch (index) {
  case 0:
    fBigBoxX = value ; 
    break ; 
  case 1:
    fBigBoxY = value ; 
    break ; 
  case 2:
    fBigBoxZ = value ; 
    break ; 
 }

}

//____________________________________________________________________________
Double_t AliPHOSvFast::SigmaE(Double_t energy)
{
  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)
{
 
  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)
{

  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) ;

  // 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
	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"
3decf5cb	35	#include "AliRun.h"
	36	#include "AliConst.h"
	37
	38	ClassImp(AliPHOSvFast)
	39
	40	//____________________________________________________________________________
	41	AliPHOSvFast::AliPHOSvFast()
	42	{
a73f33f0	43	fFastRecParticles = 0 ;
3decf5cb	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 ;
a73f33f0	65	fFastRecParticles = new FastRecParticlesList("AliPHOSFastRecParticle", 100) ;
a73f33f0	66
b1aa729d	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 ;
3decf5cb	76	}
	77
	78	//____________________________________________________________________________
	79	AliPHOSvFast::~AliPHOSvFast()
	80	{
	81
a73f33f0	82	fFastRecParticles->Delete() ;
	83	delete fFastRecParticles ;
	84	fFastRecParticles = 0 ;
3decf5cb	85
	86	}
	87
	88	//____________________________________________________________________________
a73f33f0	89	void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
	90	{
	91	new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
	92	fNRecParticles++ ;
3decf5cb	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 ) ;
3decf5cb	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());
a73f33f0	229	char *cd = strstr(opt,"R");
3decf5cb	230
a73f33f0	231	if (fFastRecParticles && gAlice->TreeR() && cd) {
a73f33f0	232	gAlice->TreeR()->Branch(branchname, &fFastRecParticles, fBufferSize);
3decf5cb	233	}
	234	}
	235
a73f33f0	236	//____________________________________________________________________________
	237	Double_t AliPHOSvFast::MakeEnergy(const Double_t energy)
	238	{
b1aa729d	239	Double_t sigma = SigmaE(energy) ;
b1aa729d	240	return fRan.Gaus(energy, sigma) ;
a73f33f0	241	}
	242
	243	//____________________________________________________________________________
b1aa729d	244	TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi)
a73f33f0	245	{
b1aa729d	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	}
a73f33f0	259
b1aa729d	260	//____________________________________________________________________________
	261	void AliPHOSvFast::MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp)
	262	{
a73f33f0	263	// get the detected type of particle
b1aa729d	264
b1aa729d	265	Int_t type = MakeType( rp ) ;
a73f33f0	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 ) ;
b1aa729d	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
a73f33f0	292	}
	293
	294	//____________________________________________________________________________
b1aa729d	295	Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp )
a73f33f0	296	{
a73f33f0	297	Int_t rv = kUNDEFINED ;
b1aa729d	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
a73f33f0	364	return rv ;
	365	}
	366
3decf5cb	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
a73f33f0	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
b1aa729d	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
3decf5cb	407	//____________________________________________________________________________
	408	void AliPHOSvFast::StepManager(void)
	409	{
	410
	411	Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
a73f33f0	412	TLorentzVector lv ;
a73f33f0	413	gMC->TrackPosition(lv) ;
b1aa729d	414	TVector3 pos = lv.Vect() ;
	415	Int_t modid ;
	416	gMC->CurrentVolID(modid);
a73f33f0	417
a73f33f0	418	// Makes a reconstructed particle from the primary particle
3decf5cb	419
a73f33f0	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
b1aa729d	426
b1aa729d	427	MakeRecParticle(modid, pos, rp) ;
3decf5cb	428
a73f33f0	429	// add the primary particle to the FastRecParticles list
b1aa729d	430
a73f33f0	431	AddRecParticle(rp) ;
	432
	433	// stop the track as soon PHOS is reached
b1aa729d	434
3decf5cb	435	gMC->StopTrack() ;
	436
	437	}
	438