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

+

+

+ Fast All Views +

+

Front View

+

+

+ Fast 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() ; }