From b1aa729da0c29bc412e01012a467c1463c469b33 Mon Sep 17 00:00:00 2001 From: schutz Date: Thu, 2 Mar 2000 11:00:21 +0000 Subject: [PATCH] Fast simulation continued: Particle identification included --- PHOS/AliPHOSRecParticle.cxx | 38 +++++----- PHOS/AliPHOSvFast.cxx | 144 ++++++++++++++++++++++++++++++++---- PHOS/AliPHOSvFast.h | 25 ++++--- 3 files changed, 166 insertions(+), 41 deletions(-) diff --git a/PHOS/AliPHOSRecParticle.cxx b/PHOS/AliPHOSRecParticle.cxx index d42a13e3ff7..1470dbc4b50 100644 --- a/PHOS/AliPHOSRecParticle.cxx +++ b/PHOS/AliPHOSRecParticle.cxx @@ -50,25 +50,25 @@ ClassImp(AliPHOSRecParticle) { fPHOSTrackSegment = new AliPHOSTrackSegment( *( rp.GetPHOSTrackSegment()) ) ; fType = rp.fType ; - fPdgCode = rp.fPdgCode; - fStatusCode = rp.fStatusCode; - fMother[0] = rp.fMother[0]; - fMother[1] = rp.fMother[1]; - fDaughter[0] = rp.fDaughter[0]; - fDaughter[1] = rp.fDaughter[1]; - fWeight = rp.fWeight; - fCalcMass = rp.fCalcMass; - fPx = rp.fPx; - fPy = rp.fPy; - fPz = rp.fPz; - fE = rp.fE; - fVx = rp.fVx; - fVy = rp.fVy; - fVz = rp.fVz; - fVt = rp.fVt; - fPolarTheta = rp.fPolarTheta; - fPolarPhi = rp.fPolarPhi; - fParticlePDG = rp.fParticlePDG; +// fPdgCode = rp.fPdgCode; +// fStatusCode = rp.fStatusCode; +// fMother[0] = rp.fMother[0]; +// fMother[1] = rp.fMother[1]; +// fDaughter[0] = rp.fDaughter[0]; +// fDaughter[1] = rp.fDaughter[1]; +// fWeight = rp.fWeight; +// fCalcMass = rp.fCalcMass; +// fPx = rp.fPx; +// fPy = rp.fPy; +// fPz = rp.fPz; +// fE = rp.fE; +// fVx = rp.fVx; +// fVy = rp.fVy; +// fVz = rp.fVz; +// fVt = rp.fVt; +// fPolarTheta = rp.fPolarTheta; +// fPolarPhi = rp.fPolarPhi; +// fParticlePDG = rp.fParticlePDG; } //____________________________________________________________________________ diff --git a/PHOS/AliPHOSvFast.cxx b/PHOS/AliPHOSvFast.cxx index d1e2b4d6274..aa8d558902d 100644 --- a/PHOS/AliPHOSvFast.cxx +++ b/PHOS/AliPHOSvFast.cxx @@ -64,10 +64,15 @@ AliPHOSvFast::AliPHOSvFast(const char *name, const char *title): fNRecParticles = 0 ; fFastRecParticles = new FastRecParticlesList("AliPHOSFastRecParticle", 100) ; - fResPara1 = 30.0 ; - fResPara2 = 0.03 ; - fResPara3 = 0.01 ; - + 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 ; } //____________________________________________________________________________ @@ -165,7 +170,6 @@ void AliPHOSvFast::CreateGeometry() 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 @@ -232,16 +236,33 @@ void AliPHOSvFast::MakeBranch(Option_t* opt) //____________________________________________________________________________ Double_t AliPHOSvFast::MakeEnergy(const Double_t energy) { - Double_t sigma = SigmaE(energy*1000.) ; - return fRan.Gaus(energy*1000., sigma)/1000. ; + Double_t sigma = SigmaE(energy) ; + return fRan.Gaus(energy, sigma) ; } //____________________________________________________________________________ -void AliPHOSvFast::MakeRecParticle(AliPHOSFastRecParticle & rp) +TVector3 AliPHOSvFast::MakePosition(const Double_t energy, const TVector3 pos, const Double_t theta, const Double_t phi) { + 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) +{ // get the detected type of particle - Int_t type = MakeType( rp.GetName() ) ; + + Int_t type = MakeType( rp ) ; rp.SetType(type) ; @@ -251,14 +272,95 @@ void AliPHOSvFast::MakeRecParticle(AliPHOSFastRecParticle & rp) rp.Momentum(momentum) ; Double_t kineticenergy = TMath::Sqrt( TMath::Power(momentum.E(), 2) - TMath::Power(rp.GetMass(), 2) ) ; Double_t modifiedkineticenergy = MakeEnergy(kineticenergy ) ; - cout << modifiedkineticenergy << endl ; - // rp.SetMomentum(tempo) ; + 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 AliPHOSvFast::MakeType(const Text_t * name) +Int_t AliPHOSvFast::MakeType(AliPHOSFastRecParticle & rp ) { Int_t rv = 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 = kGAMMA ; + else { + ran = fRan.Rndm() ; + if( ran <= 0.9498 ) + rv = kNEUTRALEM ; + else + rv = kNEUTRALHADRON ; + } + break ; + + case 2112: // it's a neutron + ran = fRan.Rndm() ; + if ( ran <= 0.9998 ) + rv = kNEUTRALHADRON ; + else + rv = kNEUTRALEM ; + break ; + + case -2112: // it's a anti-neutron + ran = fRan.Rndm() ; + if ( ran <= 0.9984 ) + rv = kNEUTRALHADRON ; + else + rv = kNEUTRALEM ; + break ; + + case 11: // it's a electron + ran = fRan.Rndm() ; + if ( ran <= 0.9996 ) + rv = kELECTRON ; + else + rv = kCHARGEDHADRON ; + break; + + case -11: // it's a electron + ran = fRan.Rndm() ; + if ( ran <= 0.9996 ) + rv = kELECTRON ; + else + rv = kCHARGEDHADRON ; + break; + + case -1: // it's a charged + ran = fRan.Rndm() ; + if ( ran <= 0.9996 ) + rv = kCHARGEDHADRON ; + else + rv = kGAMMA ; + + break ; + } + + return rv ; } @@ -292,6 +394,16 @@ Double_t AliPHOSvFast::SigmaE(Double_t energy) return rv * energy ; } +//____________________________________________________________________________ +Double_t AliPHOSvFast::SigmaP(Double_t energy, Int_t incidence) +{ + + 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) { @@ -299,6 +411,9 @@ void AliPHOSvFast::StepManager(void) Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() ); TLorentzVector lv ; gMC->TrackPosition(lv) ; + TVector3 pos = lv.Vect() ; + Int_t modid ; + gMC->CurrentVolID(modid); // Makes a reconstructed particle from the primary particle @@ -308,12 +423,15 @@ void AliPHOSvFast::StepManager(void) AliPHOSFastRecParticle rp(*part) ; // Adds the response of PHOS to the particle - MakeRecParticle(rp) ; + + MakeRecParticle(modid, pos, rp) ; // add the primary particle to the FastRecParticles list + AddRecParticle(rp) ; // stop the track as soon PHOS is reached + gMC->StopTrack() ; } diff --git a/PHOS/AliPHOSvFast.h b/PHOS/AliPHOSvFast.h index 24860250ca8..6c02ad78313 100644 --- a/PHOS/AliPHOSvFast.h +++ b/PHOS/AliPHOSvFast.h @@ -40,12 +40,15 @@ public: Int_t IsVersion(void) const { return -1 ; } void MakeBranch(Option_t* opt) ; Double_t MakeEnergy(const Double_t energy) ; // makes the detected energy - void MakeRecParticle(AliPHOSFastRecParticle & rp) ; // makes a reconstructes particle from primary - Int_t MakeType(const Text_t * name) ; // gets the detected type of particle + TVector3 MakePosition(const Double_t energy, const TVector3 pos, const Double_t th, const Double_t ph) ; + // makes the detected position + void MakeRecParticle(const Int_t modid, const TVector3 pos, AliPHOSFastRecParticle & rp) ; // makes a reconstructes particle from primary + Int_t MakeType(AliPHOSFastRecParticle & rp) ; // gets the detected type of particle FastRecParticlesList * FastRecParticles() { return fFastRecParticles ; } // gets TClonesArray of reconstructed particles - void SetBigBox(Int_t index, Float_t value) ; - Double_t SigmaE(Double_t energy) ; // calulates the energy resolution at a given Energy - virtual void StepManager(void) ; // does the tracking through PHOS and a preliminary digitalization + void SetBigBox(Int_t index, Float_t value) ; + Double_t SigmaE(Double_t energy) ; // calulates the energy resolution at a given Energy + Double_t SigmaP(Double_t energy, Int_t inc) ; // calulates the position resolution at a given Energy at a given incidence + virtual void StepManager(void) ; // does the tracking through PHOS and a preliminary digitalization private: @@ -56,10 +59,14 @@ private: AliPHOSGeometry * fGeom ; // geometry definition Int_t fNRecParticles ; // number of detected particles TRandom fRan ; // random number generator - Double_t fResPara1 ; // parameter for the energy resolution dependence ; - Double_t fResPara2 ; // parameter for the energy resolution dependence ; - Double_t fResPara3 ; // parameter for the energy resolution dependence ; - + Double_t fResPara1 ; // parameter for the energy resolution dependence + Double_t fResPara2 ; // parameter for the energy resolution dependence + Double_t fResPara3 ; // parameter for the energy resolution dependence + Double_t fPosParaA0 ; // parameter for the position resolution + Double_t fPosParaA1 ; // parameter for the position resolution + Double_t fPosParaB0 ; // parameter for the position resolution + Double_t fPosParaB1 ; // parameter for the position resolution + Double_t fPosParaB2 ; // parameter for the position resolution ClassDef(AliPHOSvFast,1) // PHOS main class , version for fast simulation -- 2.43.0