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

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

//_________________________________________________________________________
// 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 <TGeometry.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include "TClonesArray.h" 
#include <TVirtualMC.h>

// --- Standard library ---

// --- AliRoot header files ---
#include "AliPHOSFastRecParticle.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSLoader.h"
#include "AliPHOSvFast.h"
#include "AliRun.h"

ClassImp(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.)    
{
  // default ctor : initialize data member
}

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast(const char *name, const char *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)
{
  // 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 ;
  fFastRecParticles = 0 ; 

}

//____________________________________________________________________________
void AliPHOSvFast::AddRecParticle(const AliPHOSFastRecParticle & rp)
{  
  // Add a virtually reconstructed particle to the list 

  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") ;
    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 / TMath::Pi() ;
  
  for( Int_t i = 1; i <= GetGeometry()->GetNModules(); 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 = GetGeometry()->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) const
{
  // 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)
{  
  // Create new branch in the current reconstructed Root Tree
  AliDetector::MakeBranch(opt);
  const char *cd = strstr(opt,"R");
  
  if (fFastRecParticles && fLoader->TreeR() && cd) {
    MakeBranchInTree(fLoader->TreeR(), GetName(), &fFastRecParticles, fBufferSize, 0);
  }
}
//____________________________________________________________________________

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(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 ) ;
  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(Int_t modid, 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 + GetGeometry()->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() ; 

  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              // NB: ALL VALUES SHOLD BE CHECKED !!!!
    ran = fRan.Rndm() ; 
    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 =  AliPHOSFastRecParticle::kNEUTRALHASLOW ; 
    else 
      rv = AliPHOSFastRecParticle::kNEUTRALEMSLOW ; 
    break ; 
    
  case -2112: // it's a anti-neutron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9984 )
      rv =  AliPHOSFastRecParticle::kNEUTRALHASLOW ; 
    else 
      rv =  AliPHOSFastRecParticle::kNEUTRALEMSLOW ; 
    break ; 
    
  case 11:    // it's a electron
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9996 )
      rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ; 
    else 
      rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ; 
    break; 

  case -11:   // it's a positon
    ran = fRan.Rndm() ; 
    if ( ran <= 0.9996 )
      rv =  AliPHOSFastRecParticle::kCHARGEDEMFAST ; 
    else 
      rv =  AliPHOSFastRecParticle::kCHARGEDHAFAST ; 
    break; 

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

    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, Double_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 

  TLorentzVector lv ; 
  gMC->TrackPosition(lv) ;
  TVector3 pos = lv.Vect() ; 
  Int_t modid  ; 
  gMC->CurrentVolID(modid);
  
  Float_t energy = gMC->Etot() ; //Total energy of current track

  //Calculating mass of current particle
  TDatabasePDG * pdg = TDatabasePDG::Instance() ;
  TParticlePDG * partPDG = pdg->GetParticle(gMC->TrackPid()) ;
  Float_t mass = partPDG->Mass() ;

  if(energy > mass){
    pos.SetMag(TMath::Sqrt(energy*energy-mass*mass)) ;
    TLorentzVector pTrack(pos, energy) ;  

    TParticle * part = new TParticle(gMC->TrackPid(), 0,-1,-1,-1,-1, pTrack, lv)  ;
        
    AliPHOSFastRecParticle rp(*part) ;

    // Adds the response of PHOS to the particle
    MakeRecParticle(modid, pos, rp) ;
    
    // add the `track' particle to the FastRecParticles list
  
    AddRecParticle(rp) ;

    part->Delete() ;
  }
  // stop the track as soon PHOS is reached
  
  gMC->StopTrack() ; 

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