1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////
20 // Segment for reconstruction
25 // two hits for reconstruction in the two chambers of one station
27 ///////////////////////////////////////////////////////////
30 #include "AliMUONChamber.h"
31 #include "AliMUONHitForRec.h"
32 #include "AliMUONSegment.h"
33 #include "AliMUONTrackParam.h"
34 #include "AliRun.h" // for gAlice
36 ClassImp(AliMUONSegment) // Class implementation in ROOT context
38 //__________________________________________________________________________
39 AliMUONSegment::AliMUONSegment()
41 // Default constructor
42 fHitForRecPtr1 = 0; // pointer to HitForRec in first chamber
43 fHitForRecPtr2 = 0; // pointer to HitForRec in second chamber
45 fBendingCoor = 0.0; // Coordinate in bending plane
46 fBendingSlope = 0.0; // Slope in bending plane
47 // Covariance in bending plane:
48 fBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
49 fBendingSlopeReso2 = 0.0; // Covariance(slope)
50 fBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
51 fBendingImpact = 0.0; // Impact parameter in bending plane
53 fNonBendingCoor = 0.0; // Coordinate in non bending plane
54 fNonBendingSlope = 0.0; // Slope in non bending plane
55 // Covariance in non bending plane:
56 fNonBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
57 fNonBendingSlopeReso2 = 0.0; // Covariance(slope)
58 fNonBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
59 fNonBendingImpact = 0.0; // Impact parameter in non bending plane
60 fInTrack = kFALSE; // TRUE if segment belongs to one track
63 //__________________________________________________________________________
64 AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
66 // Constructor for AliMUONSegment from two HitForRec's,
67 // one, in the first chamber of the station, pointed to by "Hit1",
68 // the other one, in the second chamber of the station, pointed to by "Hit1".
69 // Fills the pointers to both hits,
70 // the slope, the covariance for (coordinate in first chamber, slope),
71 // and the impact parameter at vertex (Z=0),
72 // in bending and non bending planes.
73 // Puts the "fInTrack" flag to "kFALSE".
75 // pointers to HitForRec's
76 fHitForRecPtr1 = Hit1;
77 fHitForRecPtr2 = Hit2;
78 dz = Hit1->GetZ() - Hit2->GetZ();
80 fBendingCoor = Hit1->GetBendingCoor();
81 fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
82 fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
83 fBendingCoorReso2 = Hit1->GetBendingReso2();
84 fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
85 Hit2->GetBendingReso2() ) / dz / dz;
86 fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
88 fNonBendingCoor = Hit1->GetNonBendingCoor();
89 fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
90 fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
91 fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
92 fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
93 Hit2->GetNonBendingReso2() ) / dz / dz;
94 fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
95 // "fInTrack" flag to "kFALSE"
100 AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
102 // Dummy copy constructor
105 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
107 // Dummy assignment operator
111 //__________________________________________________________________________
112 Int_t AliMUONSegment::Compare(const TObject* Segment) const
114 // "Compare" function to sort with increasing absolute value
115 // of the "impact parameter" in bending plane.
116 // Returns -1 (0, +1) if |impact parameter| of current Segment
117 // is smaller than (equal to, larger than) |impact parameter| of Segment
118 if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
119 < TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
121 // continuous parameter, hence no need for testing equal case
125 //__________________________________________________________________________
126 Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut)
128 // Calculate the normalized Chi2 between the current Segment (this)
129 // and the Segment pointed to by "Segment",
130 // i.e. the square deviations between the coordinates and the slopes,
131 // in both the bending and the non bending plane,
132 // divided by the variance of the same quantities and by "Sigma2Cut".
133 // Returns 5 if none of the 4 quantities is OK,
134 // something smaller than or equal to 4 otherwise.
135 // Would it be more correct to use a real chi square
136 // including the non diagonal term ????
137 Double_t chi2, chi2Max, diff, normDiff;
140 // coordinate in bending plane
141 diff = this->fBendingCoor - Segment->fBendingCoor;
142 normDiff = diff * diff /
143 (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
144 if (normDiff > 1.0) return chi2Max;
145 chi2 = chi2 + normDiff;
146 // slope in bending plane
147 diff = this->fBendingSlope - Segment->fBendingSlope;
148 normDiff = diff * diff /
149 (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
150 if (normDiff > 1.0) return chi2Max;
151 chi2 = chi2 + normDiff;
152 // coordinate in non bending plane
153 diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
154 normDiff = diff * diff /
155 (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
156 if (normDiff > 1.0) return chi2Max;
157 chi2 = chi2 + normDiff;
158 // slope in non bending plane
159 diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
160 normDiff = diff * diff /
161 (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
162 if (normDiff > 1.0) return chi2Max;
163 chi2 = chi2 + normDiff;
167 //__________________________________________________________________________
168 AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor)
170 // Extrapolates linearly the current Segment (this) to station (0..) "Station".
171 // Multiple Coulomb scattering calculated from "MCSfactor"
172 // corresponding to one chamber,
173 // with one chamber for the coordinate, two chambers for the angle,
174 // due to the arrangement in stations.
175 // Valid from station(1..) 4 to 5 or vice versa.
176 // Returns the pointer to the created AliMUONSegment object
177 // corresponding to this extrapolation.
178 // The caller has the responsibility to delete this object.
179 AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
180 // dZ from first hit of current Segment to first chamber of station "Station"
181 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
183 (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
184 // Data in bending plane
186 extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
188 extrapSegment->fBendingSlope = this->fBendingSlope;
189 // covariance, including multiple Coulomb scattering over dZ due to one chamber
190 extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
191 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
192 extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
193 extrapSegment->fBendingCoorSlopeReso2 =
194 this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
195 // Data in non bending plane
197 extrapSegment->fNonBendingCoor =
198 this->fNonBendingCoor + this->fNonBendingSlope * dZ;
200 extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
201 // covariance, including multiple Coulomb scattering over dZ due to one chamber
202 extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
203 (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
204 extrapSegment->fNonBendingSlopeReso2 =
205 this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
206 extrapSegment->fNonBendingCoorSlopeReso2 =
207 this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
208 return extrapSegment;
211 //__________________________________________________________________________
212 AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor)
214 // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
215 // Multiple Coulomb scattering calculated from "MCSfactor"
216 // corresponding to one chamber.
217 // Valid from station(1..) 4 to 5 or vice versa.
218 // Returns the pointer to the created AliMUONHitForRec object
219 // corresponding to this extrapolation.
220 // The caller has the responsibility to delete this object.
221 AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
222 // dZ from first hit of current Segment to chamber
223 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
225 (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
226 // Data in bending plane
228 extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
229 // covariance, including multiple Coulomb scattering over dZ due to one chamber
230 extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
231 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
232 // Data in non bending plane
234 extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
235 this->fNonBendingSlope * dZ);
236 // covariance, including multiple Coulomb scattering over dZ due to one chamber
238 SetNonBendingReso2(this->fNonBendingCoorReso2 +
239 (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
240 return extrapHitForRec;
243 //__________________________________________________________________________
244 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
246 // Fill data members with values calculated from the array of track parameters
247 // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
248 // of the station, respectively).
249 // Multiple Coulomb scattering is taking into account with "MCSfactor"
250 // corresponding to one chamber,
251 // with one chamber for the coordinate, two chambers for the angle,
252 // due to the arrangement in stations.
253 // Resolution coming from:
254 // coordinate in closest station at "Dz1" from current "Station",
255 // slope between closest stations, with "Dz2" interval between them,
256 // interval "Dz3" between chambers of closest station,
257 // extrapolation over "Dz1" from closest station,
258 // "InverseMomentum".
259 // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
260 // are assumed to be filled
261 // with the variance on bending and non bending coordinates.
262 // The "road" is parametrized from the old reco_muon.F
263 // with 8 cm between stations.
264 AliMUONTrackParam *param0;
265 // Double_t cReso2, sReso2;
266 // parameters to define the widths of the searching roads in station 0,1,2
267 // width = p0 + p1/ (momentum)^2
268 // station number: 0 1 2
269 // static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
270 // static Double_t p1BendingCoor[3] = { 986., 821., 446. };
271 // static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
272 // static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
273 // static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
274 // static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
275 // static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
276 // static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
278 static Double_t p0BendingCoor[3] = { 6.43e-2, 6.43e-2, 6.43e-2 };
279 static Double_t p1BendingCoor[3] = { 986., 986., 986. };
280 static Double_t p0BendingSlope[3] = { 3.6e-6, 3.6e-6, 3.6e-6 };
281 static Double_t p1BendingSlope[3] = { 1.1e-2, 1.1e-2, 1.1e-2 };
282 static Double_t p0NonBendingCoor[3] = { 0.049, 0.049, 0.049 };
283 static Double_t p1NonBendingCoor[3] = { 1444., 1444., 1444. };
284 static Double_t p0NonBendingSlope[3] = { 6.8e-4, 6.8e-4, 6.8e-4 };
285 static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
286 param0 = &(TrackParam[0]);
290 // fBendingCoor = param0->GetBendingCoor(); // coordinate
291 // fBendingSlope = param0->GetBendingSlope(); // slope
292 // cReso2 = fBendingCoorReso2;
293 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
294 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
295 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
296 // // Non bending plane
297 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
298 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
299 // cReso2 = fNonBendingCoorReso2;
300 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
301 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
302 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
304 // Coordinate and slope
306 fBendingCoor = param0->GetBendingCoor(); // coordinate
307 fBendingSlope = param0->GetBendingSlope(); // slope
309 fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
310 fNonBendingSlope = param0->GetNonBendingSlope(); // slope
313 // cReso2 and sReso2 have to be subtracted here from the parametrization
314 // because they are added in the functions "NormalizedChi2WithSegment"
315 // and "NormalizedChi2WithHitForRec"
317 // cReso2 = fBendingCoorReso2;
318 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
319 fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ; // - cReso2
320 fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
322 // cReso2 = fNonBendingCoorReso2;
323 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
324 fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
325 fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
329 // OLD function, with roads automatically calculated instead from being parametrized
330 // kept because it would be a better solution,
331 // if one can really find the right values.
332 // //__________________________________________________________________________
333 // void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
335 // // Fill data members with values calculated from the array of track parameters
336 // // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
337 // // of the station, respectively).
338 // // Multiple Coulomb scattering is taking into account with "MCSfactor"
339 // // corresponding to one chamber,
340 // // with one chamber for the coordinate, two chambers for the angle,
341 // // due to the arrangement in stations.
342 // // Resolution coming from:
343 // // coordinate in closest station at "Dz1",
344 // // slope between closest stations, with "Dz2" interval between them,
345 // // extrapolation over "Dz" from closest station.
346 // // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
347 // // are assumed to be filled
348 // // with the variance on bending and non bending coordinates.
349 // AliMUONTrackParam *param0;
350 // Double_t cReso2, sReso2;
351 // param0 = &(TrackParam[0]);
353 // fBendingCoor = param0->GetBendingCoor(); // coordinate
354 // fBendingSlope = param0->GetBendingSlope(); // slope
355 // cReso2 = fBendingCoorReso2;
356 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
357 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
358 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
359 // // Non bending plane
360 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
361 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
362 // cReso2 = fNonBendingCoorReso2;
363 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
364 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
365 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;