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 Revision 1.8 2001/04/25 14:50:42 gosset
19 Corrections to violations of coding conventions
21 Revision 1.7 2001/02/08 10:34:41 gosset
22 Add a "real" default constructor.
24 Revision 1.6 2001/02/05 14:49:29 hristov
25 Compare() declared const (R.Brun)
27 Revision 1.5 2001/01/08 11:01:02 gosset
28 Modifications used for addendum to Dimuon TDR (JP Cussonneau):
29 *. MaxBendingMomentum to make both a segment and a track (default 500)
30 *. MaxChi2 per degree of freedom to make a track (default 100)
31 *. MinBendingMomentum used also to make a track
32 and not only a segment (default 3)
33 *. wider roads for track search in stations 1 to 3
34 *. extrapolation to actual Z instead of Z(chamber) in FollowTracks
36 - limits on parameters X and Y (+/-500)
37 - covariance matrices in double precision
38 - normalization of covariance matrices before inversion
39 - suppression of Minuit printouts
40 *. correction against memory leak (delete extrapHit) in FollowTracks
41 *. RMax to 10 degrees with Z(chamber) instead of fixed values;
42 RMin and Rmax cuts suppressed in NewHitForRecFromGEANT,
43 because useless with realistic geometry
45 Revision 1.4 2000/06/30 10:15:48 gosset
46 Changes to EventReconstructor...:
47 precision fit with multiple Coulomb scattering;
48 extrapolation to vertex with Branson correction in absorber (JPC)
50 Revision 1.3 2000/06/25 13:06:39 hristov
51 Inline functions moved from *.cxx to *.h files instead of forward declarations
53 Revision 1.2 2000/06/15 07:58:48 morsch
54 Code from MUON-dev joined
56 Revision 1.1.2.4 2000/06/12 10:10:21 morsch
57 Dummy copy constructor and assignment operator added
59 Revision 1.1.2.3 2000/06/09 21:01:16 morsch
60 Make includes consistent with new file structure.
62 Revision 1.1.2.2 2000/06/09 12:58:05 gosset
63 Removed comment beginnings in Log sections of .cxx files
64 Suppressed most violations of coding rules
66 Revision 1.1.2.1 2000/06/07 14:44:53 gosset
67 Addition of files for track reconstruction in C++
70 ///////////////////////////////////////////////////////////
72 // Segment for reconstruction
77 // two hits for reconstruction in the two chambers of one station
79 ///////////////////////////////////////////////////////////
82 #include "AliMUONChamber.h"
83 #include "AliMUONHitForRec.h"
84 #include "AliMUONSegment.h"
85 #include "AliMUONTrackParam.h"
86 #include "AliRun.h" // for gAlice
88 ClassImp(AliMUONSegment) // Class implementation in ROOT context
90 //__________________________________________________________________________
91 AliMUONSegment::AliMUONSegment()
93 // Default constructor
94 fHitForRecPtr1 = 0; // pointer to HitForRec in first chamber
95 fHitForRecPtr2 = 0; // pointer to HitForRec in second chamber
97 fBendingCoor = 0.0; // Coordinate in bending plane
98 fBendingSlope = 0.0; // Slope in bending plane
99 // Covariance in bending plane:
100 fBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
101 fBendingSlopeReso2 = 0.0; // Covariance(slope)
102 fBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
103 fBendingImpact = 0.0; // Impact parameter in bending plane
104 // Non Bending plane:
105 fNonBendingCoor = 0.0; // Coordinate in non bending plane
106 fNonBendingSlope = 0.0; // Slope in non bending plane
107 // Covariance in non bending plane:
108 fNonBendingCoorReso2 = 0.0; // Covariance(coordinate C1 in first chamber)
109 fNonBendingSlopeReso2 = 0.0; // Covariance(slope)
110 fNonBendingCoorSlopeReso2 = 0.0; // Covariance(C1,slope)
111 fNonBendingImpact = 0.0; // Impact parameter in non bending plane
112 fInTrack = kFALSE; // TRUE if segment belongs to one track
115 //__________________________________________________________________________
116 AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
118 // Constructor for AliMUONSegment from two HitForRec's,
119 // one, in the first chamber of the station, pointed to by "Hit1",
120 // the other one, in the second chamber of the station, pointed to by "Hit1".
121 // Fills the pointers to both hits,
122 // the slope, the covariance for (coordinate in first chamber, slope),
123 // and the impact parameter at vertex (Z=0),
124 // in bending and non bending planes.
125 // Puts the "fInTrack" flag to "kFALSE".
127 // pointers to HitForRec's
128 fHitForRecPtr1 = Hit1;
129 fHitForRecPtr2 = Hit2;
130 dz = Hit1->GetZ() - Hit2->GetZ();
132 fBendingCoor = Hit1->GetBendingCoor();
133 fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
134 fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
135 fBendingCoorReso2 = Hit1->GetBendingReso2();
136 fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
137 Hit2->GetBendingReso2() ) / dz / dz;
138 fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
140 fNonBendingCoor = Hit1->GetNonBendingCoor();
141 fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
142 fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
143 fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
144 fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
145 Hit2->GetNonBendingReso2() ) / dz / dz;
146 fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
147 // "fInTrack" flag to "kFALSE"
152 AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
154 // Dummy copy constructor
157 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
159 // Dummy assignment operator
163 //__________________________________________________________________________
164 Int_t AliMUONSegment::Compare(const TObject* Segment) const
166 // "Compare" function to sort with increasing absolute value
167 // of the "impact parameter" in bending plane.
168 // Returns -1 (0, +1) if |impact parameter| of current Segment
169 // is smaller than (equal to, larger than) |impact parameter| of Segment
170 if (TMath::Abs(((AliMUONSegment*)this)->fBendingImpact)
171 < TMath::Abs(((AliMUONSegment*)Segment)->fBendingImpact))
173 // continuous parameter, hence no need for testing equal case
177 //__________________________________________________________________________
178 Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut)
180 // Calculate the normalized Chi2 between the current Segment (this)
181 // and the Segment pointed to by "Segment",
182 // i.e. the square deviations between the coordinates and the slopes,
183 // in both the bending and the non bending plane,
184 // divided by the variance of the same quantities and by "Sigma2Cut".
185 // Returns 5 if none of the 4 quantities is OK,
186 // something smaller than or equal to 4 otherwise.
187 // Would it be more correct to use a real chi square
188 // including the non diagonal term ????
189 Double_t chi2, chi2Max, diff, normDiff;
192 // coordinate in bending plane
193 diff = this->fBendingCoor - Segment->fBendingCoor;
194 normDiff = diff * diff /
195 (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
196 if (normDiff > 1.0) return chi2Max;
197 chi2 = chi2 + normDiff;
198 // slope in bending plane
199 diff = this->fBendingSlope - Segment->fBendingSlope;
200 normDiff = diff * diff /
201 (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
202 if (normDiff > 1.0) return chi2Max;
203 chi2 = chi2 + normDiff;
204 // coordinate in non bending plane
205 diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
206 normDiff = diff * diff /
207 (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
208 if (normDiff > 1.0) return chi2Max;
209 chi2 = chi2 + normDiff;
210 // slope in non bending plane
211 diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
212 normDiff = diff * diff /
213 (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
214 if (normDiff > 1.0) return chi2Max;
215 chi2 = chi2 + normDiff;
219 //__________________________________________________________________________
220 AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor)
222 // Extrapolates linearly the current Segment (this) to station (0..) "Station".
223 // Multiple Coulomb scattering calculated from "MCSfactor"
224 // corresponding to one chamber,
225 // with one chamber for the coordinate, two chambers for the angle,
226 // due to the arrangement in stations.
227 // Valid from station(1..) 4 to 5 or vice versa.
228 // Returns the pointer to the created AliMUONSegment object
229 // corresponding to this extrapolation.
230 // The caller has the responsibility to delete this object.
231 AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
232 // dZ from first hit of current Segment to first chamber of station "Station"
233 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
235 (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
236 // Data in bending plane
238 extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
240 extrapSegment->fBendingSlope = this->fBendingSlope;
241 // covariance, including multiple Coulomb scattering over dZ due to one chamber
242 extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
243 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
244 extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
245 extrapSegment->fBendingCoorSlopeReso2 =
246 this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
247 // Data in non bending plane
249 extrapSegment->fNonBendingCoor =
250 this->fNonBendingCoor + this->fNonBendingSlope * dZ;
252 extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
253 // covariance, including multiple Coulomb scattering over dZ due to one chamber
254 extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
255 (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
256 extrapSegment->fNonBendingSlopeReso2 =
257 this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
258 extrapSegment->fNonBendingCoorSlopeReso2 =
259 this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
260 return extrapSegment;
263 //__________________________________________________________________________
264 AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor)
266 // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
267 // Multiple Coulomb scattering calculated from "MCSfactor"
268 // corresponding to one chamber.
269 // Valid from station(1..) 4 to 5 or vice versa.
270 // Returns the pointer to the created AliMUONHitForRec object
271 // corresponding to this extrapolation.
272 // The caller has the responsibility to delete this object.
273 AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
274 // dZ from first hit of current Segment to chamber
275 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
277 (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
278 // Data in bending plane
280 extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
281 // covariance, including multiple Coulomb scattering over dZ due to one chamber
282 extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
283 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
284 // Data in non bending plane
286 extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
287 this->fNonBendingSlope * dZ);
288 // covariance, including multiple Coulomb scattering over dZ due to one chamber
290 SetNonBendingReso2(this->fNonBendingCoorReso2 +
291 (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
292 return extrapHitForRec;
295 //__________________________________________________________________________
296 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
298 // Fill data members with values calculated from the array of track parameters
299 // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
300 // of the station, respectively).
301 // Multiple Coulomb scattering is taking into account with "MCSfactor"
302 // corresponding to one chamber,
303 // with one chamber for the coordinate, two chambers for the angle,
304 // due to the arrangement in stations.
305 // Resolution coming from:
306 // coordinate in closest station at "Dz1" from current "Station",
307 // slope between closest stations, with "Dz2" interval between them,
308 // interval "Dz3" between chambers of closest station,
309 // extrapolation over "Dz1" from closest station,
310 // "InverseMomentum".
311 // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
312 // are assumed to be filled
313 // with the variance on bending and non bending coordinates.
314 // The "road" is parametrized from the old reco_muon.F
315 // with 8 cm between stations.
316 AliMUONTrackParam *param0;
317 // Double_t cReso2, sReso2;
318 // parameters to define the widths of the searching roads in station 0,1,2
319 // width = p0 + p1/ (momentum)^2
320 // station number: 0 1 2
321 // static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
322 // static Double_t p1BendingCoor[3] = { 986., 821., 446. };
323 // static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
324 // static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
325 // static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
326 // static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
327 // static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
328 // static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
330 static Double_t p0BendingCoor[3] = { 6.43e-2, 6.43e-2, 6.43e-2 };
331 static Double_t p1BendingCoor[3] = { 986., 986., 986. };
332 static Double_t p0BendingSlope[3] = { 3.6e-6, 3.6e-6, 3.6e-6 };
333 static Double_t p1BendingSlope[3] = { 1.1e-2, 1.1e-2, 1.1e-2 };
334 static Double_t p0NonBendingCoor[3] = { 0.049, 0.049, 0.049 };
335 static Double_t p1NonBendingCoor[3] = { 1444., 1444., 1444. };
336 static Double_t p0NonBendingSlope[3] = { 6.8e-4, 6.8e-4, 6.8e-4 };
337 static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
338 param0 = &(TrackParam[0]);
342 // fBendingCoor = param0->GetBendingCoor(); // coordinate
343 // fBendingSlope = param0->GetBendingSlope(); // slope
344 // cReso2 = fBendingCoorReso2;
345 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
346 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
347 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
348 // // Non bending plane
349 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
350 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
351 // cReso2 = fNonBendingCoorReso2;
352 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
353 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
354 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
356 // Coordinate and slope
358 fBendingCoor = param0->GetBendingCoor(); // coordinate
359 fBendingSlope = param0->GetBendingSlope(); // slope
361 fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
362 fNonBendingSlope = param0->GetNonBendingSlope(); // slope
365 // cReso2 and sReso2 have to be subtracted here from the parametrization
366 // because they are added in the functions "NormalizedChi2WithSegment"
367 // and "NormalizedChi2WithHitForRec"
369 // cReso2 = fBendingCoorReso2;
370 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
371 fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ; // - cReso2
372 fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
374 // cReso2 = fNonBendingCoorReso2;
375 // sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
376 fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
377 fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
381 // OLD function, with roads automatically calculated instead from being parametrized
382 // kept because it would be a better solution,
383 // if one can really find the right values.
384 // //__________________________________________________________________________
385 // void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
387 // // Fill data members with values calculated from the array of track parameters
388 // // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
389 // // of the station, respectively).
390 // // Multiple Coulomb scattering is taking into account with "MCSfactor"
391 // // corresponding to one chamber,
392 // // with one chamber for the coordinate, two chambers for the angle,
393 // // due to the arrangement in stations.
394 // // Resolution coming from:
395 // // coordinate in closest station at "Dz1",
396 // // slope between closest stations, with "Dz2" interval between them,
397 // // extrapolation over "Dz" from closest station.
398 // // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
399 // // are assumed to be filled
400 // // with the variance on bending and non bending coordinates.
401 // AliMUONTrackParam *param0;
402 // Double_t cReso2, sReso2;
403 // param0 = &(TrackParam[0]);
405 // fBendingCoor = param0->GetBendingCoor(); // coordinate
406 // fBendingSlope = param0->GetBendingSlope(); // slope
407 // cReso2 = fBendingCoorReso2;
408 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
409 // fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
410 // fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
411 // // Non bending plane
412 // fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
413 // fNonBendingSlope = param0->GetNonBendingSlope(); // slope
414 // cReso2 = fNonBendingCoorReso2;
415 // sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
416 // fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
417 // fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;