[u/mrichter/AliRoot.git] / PHOS / AliPHOSv4.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"
#include "TTree.h"
#include "TGeometry.h"
#include "TFile.h"

// --- Standard library ---

#include <stdio.h>

// --- AliRoot header files ---
#include "AliPHOSFastRecParticle.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSv4.h"
#include "AliRun.h"
#include "AliConst.h"
#include "AliMC.h"

ClassImp(AliPHOSv4)

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

  // gets an instance of the geometry parameters class  
  if (strcmp(GetTitle(),"") != 0 ) 
    fGeom =  AliPHOSGeometry::GetInstance(GetTitle(), "") ; 
     
    
    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 ; 
}

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

}

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

  new( (*fFastRecParticles)[fNRecParticles] ) AliPHOSFastRecParticle(rp) ;
  fNRecParticles++ ; 
}

//____________________________________________________________________________
void AliPHOSv4::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/AliPHOSv4AllViews.gif"> 
    </CENTER></P>
    <H4> Front View </H4>
    <P>
    <CENTER>
    <IMG Align=BOTTOM ALT="Fast Front View" SRC="../images/AliPHOSv4FrontView.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 AliPHOSv4::CreateGeometry()
{
  // Create the geometry for GEANT
  
  AliPHOSv4 *phostmp = (AliPHOSv4*)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 AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::MakeBranch(Option_t* opt, const char *file)
{  
  // Create new branch in the current reconstructed Root Tree
 
  AliDetector::MakeBranch(opt,file) ;
  
  char branchname[10];
  sprintf(branchname,"%s",GetName());
  const char *cd = strstr(opt,"R");
  
  if (fFastRecParticles && gAlice->TreeR() && cd) {
    MakeBranchInTree(gAlice->TreeR(), 
                     branchname, &fFastRecParticles, fBufferSize, file);
  }
}

//____________________________________________________________________________
Double_t AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::ResetPoints()
{
  // This overloads the method in AliDetector
  
  ResetFastRecParticles() ; 
}

//___________________________________________________________________________
void AliPHOSv4::ResetFastRecParticles()
{
  // Resets the list of virtual reconstructed particles
 
  if (fFastRecParticles) 
    fFastRecParticles->Clear() ;
  fNRecParticles = 0 ; 
}

//___________________________________________________________________________
void AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::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 AliPHOSv4::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);
  
  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) ;
    rp.SetPrimary(primary) ; 

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