/*
$Log$
+Revision 1.11 2002/03/08 17:25:36 cussonno
+Update absorber energy loss and Branson corrections : simplified functions
+BransonCorrection and TotalMomentumEnergyLoss.
+
+Revision 1.10 2001/04/25 14:50:42 gosset
+Corrections to violations of coding conventions
+
+Revision 1.9 2000/10/16 15:30:40 gosset
+TotalMomentumEnergyLoss:
+correction for change in the absorber composition (JP Cussonneau)
+
Revision 1.8 2000/10/02 21:28:09 fca
Removal of useless dependecies via forward declarations
Addition of files for track reconstruction in C++
*/
-//__________________________________________________________________________
+///////////////////////////////////////////////////
+//
+// Track parameters
+// in
+// ALICE
+// dimuon
+// spectrometer
//
-// Track parameters in ALICE dimuon spectrometer
-//__________________________________________________________________________
+///////////////////////////////////////////////////
#include <iostream.h>
#include "AliCallf77.h"
-#include "AliMUON.h"
-#include "AliMUONHitForRec.h"
-#include "AliMUONSegment.h"
+#include "AliMUON.h"
#include "AliMUONTrackParam.h"
-#include "AliMUONChamber.h"
+#include "AliMUONChamber.h"
#include "AliRun.h"
#include "AliMagF.h"
}
}
-// Inline functions for Get and Set: inline removed because it does not work !!!!
-Double_t AliMUONTrackParam::GetInverseBendingMomentum(void) {
- // Get fInverseBendingMomentum
- return fInverseBendingMomentum;}
-void AliMUONTrackParam::SetInverseBendingMomentum(Double_t InverseBendingMomentum) {
- // Set fInverseBendingMomentum
- fInverseBendingMomentum = InverseBendingMomentum;}
-Double_t AliMUONTrackParam::GetBendingSlope(void) {
- // Get fBendingSlope
- return fBendingSlope;}
-void AliMUONTrackParam::SetBendingSlope(Double_t BendingSlope) {
- // Set fBendingSlope
- fBendingSlope = BendingSlope;}
-Double_t AliMUONTrackParam::GetNonBendingSlope(void) {
- // Get fNonBendingSlope
- return fNonBendingSlope;}
-void AliMUONTrackParam::SetNonBendingSlope(Double_t NonBendingSlope) {
- // Set fNonBendingSlope
- fNonBendingSlope = NonBendingSlope;}
-Double_t AliMUONTrackParam::GetZ(void) {
- // Get fZ
- return fZ;}
-void AliMUONTrackParam::SetZ(Double_t Z) {
- // Set fZ
- fZ = Z;}
-Double_t AliMUONTrackParam::GetBendingCoor(void) {
- // Get fBendingCoor
- return fBendingCoor;}
-void AliMUONTrackParam::SetBendingCoor(Double_t BendingCoor) {
- // Set fBendingCoor
- fBendingCoor = BendingCoor;}
-Double_t AliMUONTrackParam::GetNonBendingCoor(void) {
- // Get fNonBendingCoor
- return fNonBendingCoor;}
-void AliMUONTrackParam::SetNonBendingCoor(Double_t NonBendingCoor) {
- // Set fNonBendingCoor
- fNonBendingCoor = NonBendingCoor;}
-
//__________________________________________________________________________
void AliMUONTrackParam::ExtrapToZ(Double_t Z)
{
Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
// Extrapolates track parameters upstream to the "Z" end of the front absorber
- ExtrapToZ(zAbsorber);
+ ExtrapToZ(zAbsorber); // !!!
// Makes Branson correction (multiple scattering + energy loss)
BransonCorrection();
+ // Makes a simple magnetic field correction through the absorber
+ FieldCorrection(zAbsorber);
}
+
+// Keep this version for future developments
//__________________________________________________________________________
+// void AliMUONTrackParam::BransonCorrection()
+// {
+// // Branson correction of track parameters
+// // the entry parameters have to be calculated at the end of the absorber
+// Double_t zEndAbsorber, zBP, xBP, yBP;
+// Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
+// Int_t sign;
+// // Would it be possible to calculate all that from Geant configuration ????
+// // and to get the Branson parameters from a function in ABSO module ????
+// // with an eventual contribution from other detectors like START ????
+// // Radiation lengths outer part theta > 3 degres
+// static Double_t x01[9] = { 18.8, // C (cm)
+// 10.397, // Concrete (cm)
+// 0.56, // Plomb (cm)
+// 47.26, // Polyethylene (cm)
+// 0.56, // Plomb (cm)
+// 47.26, // Polyethylene (cm)
+// 0.56, // Plomb (cm)
+// 47.26, // Polyethylene (cm)
+// 0.56 }; // Plomb (cm)
+// // inner part theta < 3 degres
+// static Double_t x02[3] = { 18.8, // C (cm)
+// 10.397, // Concrete (cm)
+// 0.35 }; // W (cm)
+// // z positions of the materials inside the absober outer part theta > 3 degres
+// static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
+// // inner part theta < 3 degres
+// static Double_t z2[4] = { 90, 315, 467, 503 };
+// static Bool_t first = kTRUE;
+// static Double_t zBP1, zBP2, rLimit;
+// // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
+// if (first) {
+// first = kFALSE;
+// Double_t aNBP = 0.0;
+// Double_t aDBP = 0.0;
+// Int_t iBound;
+
+// for (iBound = 0; iBound < 9; iBound++) {
+// aNBP = aNBP +
+// (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
+// z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound];
+// aDBP = aDBP +
+// (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound];
+// }
+// zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
+// aNBP = 0.0;
+// aDBP = 0.0;
+// for (iBound = 0; iBound < 3; iBound++) {
+// aNBP = aNBP +
+// (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
+// z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound];
+// aDBP = aDBP +
+// (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
+// }
+// zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
+// rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
+// }
+
+// pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+// sign = 1;
+// if (fInverseBendingMomentum < 0) sign = -1;
+// pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
+// pX = pZ * fNonBendingSlope;
+// pY = pZ * fBendingSlope;
+// pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
+// xEndAbsorber = fNonBendingCoor;
+// yEndAbsorber = fBendingCoor;
+// radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
+
+// if (radiusEndAbsorber2 > rLimit*rLimit) {
+// zEndAbsorber = z1[9];
+// zBP = zBP1;
+// } else {
+// zEndAbsorber = z2[3];
+// zBP = zBP2;
+// }
+
+// xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
+// yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
+
+// // new parameters after Branson and energy loss corrections
+// pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
+// pX = pZ * xBP / zBP;
+// pY = pZ * yBP / zBP;
+// fBendingSlope = pY / pZ;
+// fNonBendingSlope = pX / pZ;
+
+// pT = TMath::Sqrt(pX * pX + pY * pY);
+// theta = TMath::ATan2(pT, pZ);
+// pTotal =
+// TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
+
+// fInverseBendingMomentum = (sign / pTotal) *
+// TMath::Sqrt(1.0 +
+// fBendingSlope * fBendingSlope +
+// fNonBendingSlope * fNonBendingSlope) /
+// TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
+
+// // vertex position at (0,0,0)
+// // should be taken from vertex measurement ???
+// fBendingCoor = 0.0;
+// fNonBendingCoor = 0;
+// fZ= 0;
+// }
+
void AliMUONTrackParam::BransonCorrection()
{
// Branson correction of track parameters
// the entry parameters have to be calculated at the end of the absorber
- Double_t zEndAbsorber, zBP, xBP, yBP;
- Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
- Int_t sign;
+ // simplified version: the z positions of Branson's planes are no longer calculated
+ // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C
+ // to test this correction.
// Would it be possible to calculate all that from Geant configuration ????
// and to get the Branson parameters from a function in ABSO module ????
// with an eventual contribution from other detectors like START ????
- // Radiation lengths outer part theta > 3 degres
- static Double_t x01[9] = { 18.8, // C (cm)
- 10.397, // Concrete (cm)
- 0.56, // Plomb (cm)
- 47.26, // Polyethylene (cm)
- 0.56, // Plomb (cm)
- 47.26, // Polyethylene (cm)
- 0.56, // Plomb (cm)
- 47.26, // Polyethylene (cm)
- 0.56 }; // Plomb (cm)
- // inner part theta < 3 degres
- static Double_t x02[3] = { 18.8, // C (cm)
- 10.397, // Concrete (cm)
- 0.35 }; // W (cm)
- // z positions of the materials inside the absober outer part theta > 3 degres
- static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
- // inner part theta < 3 degres
- static Double_t z2[4] = { 90, 315, 467, 503 };
+ Double_t zBP, xBP, yBP;
+ Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
+ Int_t sign;
static Bool_t first = kTRUE;
- static Double_t zBP1, zBP2, rLimit;
- // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
+ static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber;
+ // zBP1 for outer part and zBP2 for inner part (only at the first call)
if (first) {
first = kFALSE;
- Double_t aNBP = 0.0;
- Double_t aDBP = 0.0;
- Int_t iBound;
-
- for (iBound = 0; iBound < 9; iBound++) {
- aNBP = aNBP +
- (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
- z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound];
- aDBP = aDBP +
- (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound];
- }
- zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
- aNBP = 0.0;
- aDBP = 0.0;
- for (iBound = 0; iBound < 3; iBound++) {
- aNBP = aNBP +
- (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
- z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound];
- aDBP = aDBP +
- (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
- }
- zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
- rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
+
+ zEndAbsorber = 503;
+ thetaLimit = 3.0 * (TMath::Pi()) / 180.;
+ rLimit = zEndAbsorber * TMath::Tan(thetaLimit);
+ zBP1 = 450; // values close to those calculated with EvalAbso.C
+ zBP2 = 480;
}
pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
if (radiusEndAbsorber2 > rLimit*rLimit) {
- zEndAbsorber = z1[9];
zBP = zBP1;
} else {
- zEndAbsorber = z2[3];
zBP = zBP2;
}
yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
// new parameters after Branson and energy loss corrections
- pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
- pX = pZ * xBP / zBP;
- pY = pZ * yBP / zBP;
+// Float_t zSmear = zBP - gRandom->Gaus(0.,2.); // !!! possible smearing of Z vertex position
+ Float_t zSmear = zBP;
+
+ pZ = pTotal * zSmear / TMath::Sqrt(xBP * xBP + yBP * yBP + zSmear * zSmear);
+ pX = pZ * xBP / zSmear;
+ pY = pZ * yBP / zSmear;
fBendingSlope = pY / pZ;
fNonBendingSlope = pX / pZ;
pT = TMath::Sqrt(pX * pX + pY * pY);
theta = TMath::ATan2(pT, pZ);
- pTotal =
- TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
+ pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta);
fInverseBendingMomentum = (sign / pTotal) *
TMath::Sqrt(1.0 +
}
//__________________________________________________________________________
-Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
+Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta)
{
// Returns the total momentum corrected from energy loss in the front absorber
+ // One can use the macros MUONTestAbso.C and DrawTestAbso.C
+ // to test this correction.
+ // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002)
Double_t deltaP, pTotalCorrected;
- Double_t radiusEndAbsorber2 =
- xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
- // Parametrization to be redone according to change of absorber material ????
+ // Parametrization to be redone according to change of absorber material ????
// See remark in function BransonCorrection !!!!
// The name is not so good, and there are many arguments !!!!
- if (radiusEndAbsorber2 < rLimit * rLimit) {
- if (pTotal < 15) {
- deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
+ if (theta < thetaLimit ) {
+ if (pTotal < 20) {
+ deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal;
} else {
- deltaP = 3.0643 + 0.01346 *pTotal;
+ deltaP = 3.0714 + 0.011767 *pTotal;
}
} else {
- if (pTotal < 15) {
- deltaP = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
+ if (pTotal < 20) {
+ deltaP = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal;
} else {
- deltaP = 2.407 + 0.00702 * pTotal;
+ deltaP = 2.6069 + 0.0051705 * pTotal;
}
}
- deltaP = 0.88 * deltaP; // !!!! changes in the absorber composition ????
pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
return pTotalCorrected;
}
+ //__________________________________________________________________________
+void AliMUONTrackParam::FieldCorrection(Double_t Z)
+{
+ //
+ // Correction of the effect of the magnetic field in the absorber
+ // Assume a constant field along Z axis.
+
+ Float_t b[3],x[3];
+ Double_t bZ;
+ Double_t pYZ,pX,pY,pZ,pT;
+ Double_t pXNew,pYNew;
+ Double_t c;
+
+ pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+ c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge
+
+ pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
+ pX = pZ * fNonBendingSlope;
+ pY = pZ * fBendingSlope;
+ pT = TMath::Sqrt(pX*pX+pY*pY);
+
+ if (pZ <= 0) return;
+ x[2] = Z/2;
+ x[0] = x[2]*fNonBendingSlope;
+ x[1] = x[2]*fBendingSlope;
+
+ // Take magn. field value at position x.
+ gAlice->Field()->Field(x, b);
+ bZ = b[2];
+
+ // Transverse momentum rotation
+ // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ.
+ Double_t phiShift = c*0.436*0.0003*bZ*Z/pZ;
+
+ // Rotate momentum around Z axis.
+ pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift);
+ pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift);
+
+ fBendingSlope = pYNew / pZ;
+ fNonBendingSlope = pXNew / pZ;
+
+ fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ);
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff