X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSvFast.cxx;h=53b0c2332e87997b91f0e22fc957d663fe9f0fbc;hb=8dacd1bb7c15ea67372df99783ae9dd8cef668f8;hp=08c5216c68487ee4b15610ecd102000b226e5587;hpb=3decf5cbd9be6a136c9b4d5ca807a15be3c09a97;p=u%2Fmrichter%2FAliRoot.git
diff --git a/PHOS/AliPHOSvFast.cxx b/PHOS/AliPHOSvFast.cxx
index 08c5216c684..53b0c2332e8 100644
--- a/PHOS/AliPHOSvFast.cxx
+++ b/PHOS/AliPHOSvFast.cxx
@@ -1,4 +1,4 @@
-/**************************************************************************
+ /**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
@@ -13,85 +13,137 @@
* 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.
+ *
+ */
+
//_________________________________________________________________________
-// Manager class for PHOS version for fast simulations
-//*-- Author : Y. Schutz SUBATECH
-//////////////////////////////////////////////////////////////////////////////
+// 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
+#include
+#include
+#include
+#include "TClonesArray.h"
+#include
// --- Standard library ---
-#include
-#include
-
// --- AliRoot header files ---
-
+#include "AliPHOSFastRecParticle.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSLoader.h"
#include "AliPHOSvFast.h"
-#include "AliPHOSReconstructioner.h"
#include "AliRun.h"
-#include "AliConst.h"
ClassImp(AliPHOSvFast)
-//____________________________________________________________________________
-AliPHOSvFast::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.)
{
- fRecParticles = 0 ;
- fNRecParticles = 0 ;
+ // default ctor : initialize data member
}
//____________________________________________________________________________
AliPHOSvFast::AliPHOSvFast(const char *name, const char *title):
- AliPHOS(name,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)
{
- // gets an instance of the geometry parameters class
-
- fGeom = AliPHOSGeometry::GetInstance(title, "") ;
-
- if (fGeom->IsInitialized() )
- cout << "AliPHOSvFast : PHOS geometry intialized for " << fGeom->GetName() << endl ;
- else
- cout << "AliPHOSvFast : PHOS geometry initialization failed !" << endl ;
-
- SetBigBox(0, fGeom->GetOuterBoxSize(0) ) ;
- SetBigBox(1, fGeom->GetOuterBoxSize(1) + fGeom->GetPPSDBoxSize(1) ) ;
- SetBigBox(2, fGeom->GetOuterBoxSize(0) );
-
- fNRecParticles = 0 ;
+ // 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
- fRecParticles->Delete() ;
- delete fRecParticles ;
- fRecParticles = 0 ;
+ fFastRecParticles->Delete() ;
+ delete fFastRecParticles ;
+ fFastRecParticles = 0 ;
}
//____________________________________________________________________________
-void AliPHOSvFast::AddRecParticle(Int_t primary)
-{
- TClonesArray * particlelist = gAlice->Particles() ;
- TParticle * part = (TParticle *)particlelist->At(primary) ;
- cout << " AliPHOSvFast::AddRecParticle " << part->GetName() << endl ;
+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
+ /*
+
+ PHOS FAST in ALICE displayed by root
+
+ All Views
+
+
+
+
+ Front View
+
+
+
+
+ */
+ //END_HTML
+
const Int_t kColorPHOS = kRed ;
- Double_t const kRADDEG = 180.0 / kPI ;
+ Double_t const kRADDEG = 180.0 / TMath::Pi() ;
new TBRIK( "BigBox", "PHOS box", "void", GetBigBox(0)/2,
GetBigBox(1)/2,
@@ -99,17 +151,17 @@ void AliPHOSvFast::BuildGeometry()
// position PHOS into ALICE
- Float_t r = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.0 ;
+ Float_t r = GetGeometry()->GetIPtoCrystalSurface() + GetBigBox(1) / 2.0 ;
Int_t number = 988 ;
- Float_t pphi = TMath::ATan( GetBigBox(0) / ( 2.0 * fGeom->GetIPtoOuterCoverDistance() ) ) ;
+ 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 <= fGeom->GetNModules(); i++ ) {
- Float_t angle = pphi * 2 * ( i - fGeom->GetNModules() / 2.0 - 0.5 ) ;
+ 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();
@@ -127,9 +179,10 @@ void AliPHOSvFast::BuildGeometry()
//____________________________________________________________________________
void AliPHOSvFast::CreateGeometry()
{
-
+ // Create the geometry for GEANT
+
AliPHOSvFast *phostmp = (AliPHOSvFast*)gAlice->GetModule("PHOS") ;
-
+
if ( phostmp == NULL ) {
fprintf(stderr, "PHOS detector not found!\n") ;
@@ -150,18 +203,17 @@ void AliPHOSvFast::CreateGeometry()
// --- Position PHOS mdules in ALICE setup ---
Int_t idrotm[99] ;
- Double_t const kRADDEG = 180.0 / kPI ;
+ Double_t const kRADDEG = 180.0 / TMath::Pi() ;
- for( Int_t i = 1; i <= fGeom->GetNModules(); i++ ) {
+ for( Int_t i = 1; i <= GetGeometry()->GetNModules(); i++ ) {
- Float_t angle = fGeom->GetPHOSAngle(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 = fGeom->GetIPtoOuterCoverDistance() + GetBigBox(1) / 2.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
@@ -172,7 +224,8 @@ void AliPHOSvFast::CreateGeometry()
//____________________________________________________________________________
void AliPHOSvFast::Init(void)
{
-
+ // Prints out an information message
+
Int_t i;
printf("\n");
@@ -189,8 +242,10 @@ void AliPHOSvFast::Init(void)
}
//___________________________________________________________________________
-Float_t AliPHOSvFast::GetBigBox(Int_t index)
+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) {
@@ -206,30 +261,180 @@ Float_t AliPHOSvFast::GetBigBox(Int_t index)
}
return rv ;
}
-
//___________________________________________________________________________
+
void AliPHOSvFast::MakeBranch(Option_t* opt)
{
- //
- // Create a new branch in the current Root Tree
- // The branch of fHits is automatically split
- //
- AliDetector::MakeBranch(opt) ;
+ // 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() ;
- char branchname[10];
- sprintf(branchname,"%s",GetName());
- char *cd = strstr(opt,"D");
+ 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
- if (fDigits && gAlice->TreeD() && cd) {
- gAlice->TreeD()->Branch(branchname, &fRecParticles, fBufferSize);
- // printf("* AliPHOS::MakeBranch * Making Branch %s for RecParticles \n",branchname);
+ 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::SetBigBox(Int_t index, Float_t value)
+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 ;
@@ -244,17 +449,68 @@ void AliPHOSvFast::SetBigBox(Int_t index, Float_t value)
}
+//____________________________________________________________________________
+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() ;
- Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
- TString name = fGeom->GetName() ;
+ if(energy > mass){
+ pos.SetMag(TMath::Sqrt(energy*energy-mass*mass)) ;
+ TLorentzVector pTrack(pos, energy) ;
- // add the primary particle to the RecParticles list
+ 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
- AddRecParticle(primary);
- fNRecParticles++ ;
gMC->StopTrack() ;
}