[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"
#include "TGeometry.h"

// --- Standard library ---

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

ClassImp(AliPHOSvFast)

//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast() : AliPHOS()
{
  // default ctor : initialize data member
   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. ;    

}

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

  
  // create the Getter 
  AliPHOSGetter::GetInstance(gDirectory->GetName(), 0) ; 
  
  SetBigBox(0, GetGeometry()->GetOuterBoxSize(0) ) ;
  SetBigBox(1, GetGeometry()->GetOuterBoxSize(3) + GetGeometry()->GetCPVBoxSize(1) ) ; 
  SetBigBox(2, GetGeometry()->GetOuterBoxSize(2) ); 
  
  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(const AliPHOSvFast & pv):
  AliPHOS(pv.GetName(), pv.GetTitle())
{
  // cpy ctor

  
  // create the Getter 
  AliPHOSGetter::GetInstance(gDirectory->GetName(), 0) ; 
  
  SetBigBox(0, pv.fBigBoxX) ;
  SetBigBox(1, pv.fBigBoxY) ; 
  SetBigBox(2, pv.fBigBoxZ) ; 
  
  fNRecParticles = pv.fNRecParticles; 
  fFastRecParticles = new TClonesArray( *(pv.fFastRecParticles) ) ;
  
  fResPara1 = pv.fResPara1 ;    
  fResPara2 = pv.fResPara2 ; 
  fResPara3 = pv.fResPara3 ; 
  
  fPosParaA0 = pv.fPosParaA0 ;    
  fPosParaA1 = pv.fPosParaA1 ;  
  fPosParaB0 = pv.fPosParaB0 ;   
  fPosParaB1 = pv.fPosParaB1 ; 
  fPosParaB2 = pv.fPosParaB2 ; 
}

//____________________________________________________________________________
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 = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
  Int_t number = 988 ; 
  Float_t pphi =  TMath::ATan( GetBigBox(0)  / ( 2.0 * GetGeometry()->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 <= GetGeometry()->GetNModules(); i++ ) { 
   Float_t angle = pphi * 2 * ( i - GetGeometry()->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 <= 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, 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 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 + 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 ") ;

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