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 ///////////////////////////////////////////////////////////
29 #include "AliMUONSegment.h"
31 #include "AliMUONChamber.h"
32 #include "AliMUONHitForRec.h"
33 #include "AliMUONTrackParam.h"
34 #include "AliRun.h" // for gAlice
37 ClassImp(AliMUONSegment) // Class implementation in ROOT context
39 //__________________________________________________________________________
40 AliMUONSegment::AliMUONSegment()
43 // Default constructor
44 fHitForRecPtr1 = 0; // pointer to HitForRec in first chamber
45 fHitForRecPtr2 = 0; // pointer to HitForRec in second chamber
47 fBendingCoor = 0.0; // Coordinate in bending plane
48 fBendingSlope = 0.0; // Slope in bending plane
49 // Covariance in bending plane:
50 fBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
51 fBendingSlopeReso2 = 0.0; // Covariance(slope)
52 fBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
53 fBendingImpact = 0.0; // Impact parameter in bending plane
55 fNonBendingCoor = 0.0; // Coordinate in non bending plane
56 fNonBendingSlope = 0.0; // Slope in non bending plane
57 // Covariance in non bending plane:
58 fNonBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
59 fNonBendingSlopeReso2 = 0.0; // Covariance(slope)
60 fNonBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
61 fZ = 0.0; // z in first plane
62 fNonBendingImpact = 0.0; // Impact parameter in non bending plane
63 fInTrack = kFALSE; // TRUE if segment belongs to one track
66 //__________________________________________________________________________
67 AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
70 // Constructor for AliMUONSegment from two HitForRec's,
71 // one, in the first chamber of the station, pointed to by "Hit1",
72 // the other one, in the second chamber of the station, pointed to by "Hit1".
73 // Fills the pointers to both hits,
74 // the slope, the covariance for (coordinate in first chamber, slope),
75 // and the impact parameter at vertex (Z=0),
76 // in bending and non bending planes.
77 // Puts the "fInTrack" flag to "kFALSE".
79 // pointers to HitForRec's
80 fHitForRecPtr1 = Hit1;
81 fHitForRecPtr2 = Hit2;
82 dz = Hit1->GetZ() - Hit2->GetZ();
85 fBendingCoor = Hit1->GetBendingCoor();
86 fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
87 fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
88 fBendingCoorReso2 = Hit1->GetBendingReso2();
89 fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
90 Hit2->GetBendingReso2() ) / dz / dz;
91 fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
93 fNonBendingCoor = Hit1->GetNonBendingCoor();
94 fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
95 fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
96 fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
97 fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
98 Hit2->GetNonBendingReso2() ) / dz / dz;
99 fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
100 // "fInTrack" flag to "kFALSE"
105 AliMUONSegment::AliMUONSegment (const AliMUONSegment& theMUONSegment)
106 : TObject(theMUONSegment)
108 // Protected copy constructor
110 AliFatal("Not implemented.");
113 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& rhs)
115 // Protected assignement operator
117 if (this == &rhs) return *this;
119 AliFatal("Not implemented.");
124 //__________________________________________________________________________
125 Int_t AliMUONSegment::Compare(const TObject* Segment) const
127 // "Compare" function to sort with increasing absolute value
128 // of the "impact parameter" in bending plane.
129 // Returns -1 (0, +1) if |impact parameter| of current Segment
130 // is smaller than (equal to, larger than) |impact parameter| of Segment
131 if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
132 < TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
134 // continuous parameter, hence no need for testing equal case
138 //__________________________________________________________________________
139 Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut) const
141 // Calculate the normalized Chi2 between the current Segment (this)
142 // and the Segment pointed to by "Segment",
143 // i.e. the square deviations between the coordinates and the slopes,
144 // in both the bending and the non bending plane,
145 // divided by the variance of the same quantities and by "Sigma2Cut".
146 // Returns 5 if none of the 4 quantities is OK,
147 // something smaller than or equal to 4 otherwise.
148 // Would it be more correct to use a real chi square
149 // including the non diagonal term ????
150 Double_t chi2, chi2Max, diff, normDiff;
153 // coordinate in bending plane
154 diff = this->fBendingCoor - Segment->fBendingCoor;
155 normDiff = diff * diff /
156 (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
157 if (normDiff > 1.0) return chi2Max;
158 chi2 = chi2 + normDiff;
159 // slope in bending plane
160 diff = this->fBendingSlope - Segment->fBendingSlope;
161 normDiff = diff * diff /
162 (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
163 if (normDiff > 1.0) return chi2Max;
164 chi2 = chi2 + normDiff;
165 // coordinate in non bending plane
166 diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
167 normDiff = diff * diff /
168 (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
169 if (normDiff > 1.0) return chi2Max;
170 chi2 = chi2 + normDiff;
171 // slope in non bending plane
172 diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
173 normDiff = diff * diff /
174 (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
175 if (normDiff > 1.0) return chi2Max;
176 chi2 = chi2 + normDiff;
180 //__________________________________________________________________________
181 AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation ( Double_t z, Double_t MCSfactor) const
183 // Extrapolates linearly the current Segment (this) to station (0..) "Station".
184 // Multiple Coulomb scattering calculated from "MCSfactor"
185 // corresponding to one chamber,
186 // with one chamber for the coordinate, two chambers for the angle,
187 // due to the arrangement in stations.
188 // Valid from station(1..) 4 to 5 or vice versa.
189 // Returns the pointer to the created AliMUONSegment object
190 // corresponding to this extrapolation.
191 // The caller has the responsibility to delete this object.
192 AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
193 // dZ from first hit of current Segment to first chamber of station "Station"
194 Double_t dZ = z - this->GetZ();
195 // Data in bending plane
196 extrapSegment->fZ = z;
198 extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
200 extrapSegment->fBendingSlope = this->fBendingSlope;
201 // covariance, including multiple Coulomb scattering over dZ due to one chamber
202 extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
203 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
204 extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
205 extrapSegment->fBendingCoorSlopeReso2 =
206 this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
207 // Data in non bending plane
209 extrapSegment->fNonBendingCoor =
210 this->fNonBendingCoor + this->fNonBendingSlope * dZ;
212 extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
213 // covariance, including multiple Coulomb scattering over dZ due to one chamber
214 extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
215 (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
216 extrapSegment->fNonBendingSlopeReso2 =
217 this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
218 extrapSegment->fNonBendingCoorSlopeReso2 =
219 this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
220 return extrapSegment;
223 //__________________________________________________________________________
224 AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber ( Double_t z, Double_t MCSfactor) const
226 // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
227 // Multiple Coulomb scattering calculated from "MCSfactor"
228 // corresponding to one chamber.
229 // Valid from station(1..) 4 to 5 or vice versa.
230 // Returns the pointer to the created AliMUONHitForRec object
231 // corresponding to this extrapolation.
232 // The caller has the responsibility to delete this object.
233 AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
234 // dZ from first hit of current Segment to chamber
235 Double_t dZ = z - this->GetZ();
236 // Data in bending plane
237 extrapHitForRec->SetZ(z);
239 extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
240 // covariance, including multiple Coulomb scattering over dZ due to one chamber
241 extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
242 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
243 // Data in non bending plane
245 extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
246 this->fNonBendingSlope * dZ);
247 // covariance, including multiple Coulomb scattering over dZ due to one chamber
249 SetNonBendingReso2(this->fNonBendingCoorReso2 +
250 (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
251 return extrapHitForRec;
254 //__________________________________________________________________________
255 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t /*MCSfactor*/, Double_t /*Dz1*/, Double_t /*Dz2*/, Double_t /*Dz3*/, Int_t Station, Double_t InverseMomentum)
257 // Fill data members with values calculated from the array of track parameters
258 // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
259 // of the station, respectively).
260 // Multiple Coulomb scattering is taking into account with "MCSfactor"
261 // corresponding to one chamber,
262 // with one chamber for the coordinate, two chambers for the angle,
263 // due to the arrangement in stations.
264 // Resolution coming from:
265 // coordinate in closest station at "Dz1" from current "Station",
266 // slope between closest stations, with "Dz2" interval between them,
267 // interval "Dz3" between chambers of closest station,
268 // extrapolation over "Dz1" from closest station,
269 // "InverseMomentum".
270 // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
271 // are assumed to be filled
272 // with the variance on bending and non bending coordinates.
273 // The "road" is parametrized from the old reco_muon.F
274 // with 8 cm between stations.
275 AliMUONTrackParam *param0;
276 // Double_t cReso2, sReso2;
277 // parameters to define the widths of the searching roads in station 0,1,2
278 // width = p0 + p1/ (momentum)^2
279 // station number: 0 1 2
280 // static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
281 // static Double_t p1BendingCoor[3] = { 986., 821., 446. };
282 // static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
283 // static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
284 // static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
285 // static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
286 // static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
287 // static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
289 static Double_t p0BendingCoor[3] = { 6.43e-2, 6.43e-2, 6.43e-2 };
290 static Double_t p1BendingCoor[3] = { 986., 986., 986. };
291 static Double_t p0BendingSlope[3] = { 3.6e-6, 3.6e-6, 3.6e-6 };
292 static Double_t p1BendingSlope[3] = { 1.1e-2, 1.1e-2, 1.1e-2 };
293 static Double_t p0NonBendingCoor[3] = { 0.049, 0.049, 0.049 };
294 static Double_t p1NonBendingCoor[3] = { 1444., 1444., 1444. };
295 static Double_t p0NonBendingSlope[3] = { 6.8e-4, 6.8e-4, 6.8e-4 };
296 static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
297 param0 = &(TrackParam[0]);
301 // fBendingCoor = param0->GetBendingCoor(); // coordinate
302 // fBendingSlope = param0->GetBendingSlope(); // slope
303 // cReso2 = fBendingCoorReso2;
304 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
305 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
306 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
307 // // Non bending plane
308 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
309 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
310 // cReso2 = fNonBendingCoorReso2;
311 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
312 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
313 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
315 // Coordinate and slope
317 fBendingCoor = param0->GetBendingCoor(); // coordinate
318 fBendingSlope = param0->GetBendingSlope(); // slope
320 fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
321 fNonBendingSlope = param0->GetNonBendingSlope(); // slope
323 fZ = param0->GetZ(); // z
326 // cReso2 and sReso2 have to be subtracted here from the parametrization
327 // because they are added in the functions "NormalizedChi2WithSegment"
328 // and "NormalizedChi2WithHitForRec"
330 // cReso2 = fBendingCoorReso2;
331 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
332 fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ; // - cReso2
333 fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
335 // cReso2 = fNonBendingCoorReso2;
336 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
337 fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
338 fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
342 // OLD function, with roads automatically calculated instead from being parametrized
343 // kept because it would be a better solution,
344 // if one can really find the right values.
345 // //__________________________________________________________________________
346 // void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
348 // // Fill data members with values calculated from the array of track parameters
349 // // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
350 // // of the station, respectively).
351 // // Multiple Coulomb scattering is taking into account with "MCSfactor"
352 // // corresponding to one chamber,
353 // // with one chamber for the coordinate, two chambers for the angle,
354 // // due to the arrangement in stations.
355 // // Resolution coming from:
356 // // coordinate in closest station at "Dz1",
357 // // slope between closest stations, with "Dz2" interval between them,
358 // // extrapolation over "Dz" from closest station.
359 // // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
360 // // are assumed to be filled
361 // // with the variance on bending and non bending coordinates.
362 // AliMUONTrackParam *param0;
363 // Double_t cReso2, sReso2;
364 // param0 = &(TrackParam[0]);
366 // fBendingCoor = param0->GetBendingCoor(); // coordinate
367 // fBendingSlope = param0->GetBendingSlope(); // slope
368 // cReso2 = fBendingCoorReso2;
369 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
370 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
371 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
372 // // Non bending plane
373 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
374 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
375 // cReso2 = fNonBendingCoorReso2;
376 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
377 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
378 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;