Changes to EventReconstructor...:
[u/mrichter/AliRoot.git] / MUON / AliMUONSegment.cxx
CommitLineData
a9e2aefa 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/*
17$Log$
04b5ea16 18Revision 1.3 2000/06/25 13:06:39 hristov
19Inline functions moved from *.cxx to *.h files instead of forward declarations
20
9b03f36b 21Revision 1.2 2000/06/15 07:58:48 morsch
22Code from MUON-dev joined
23
a9e2aefa 24Revision 1.1.2.4 2000/06/12 10:10:21 morsch
25Dummy copy constructor and assignment operator added
26
27Revision 1.1.2.3 2000/06/09 21:01:16 morsch
28Make includes consistent with new file structure.
29
30Revision 1.1.2.2 2000/06/09 12:58:05 gosset
31Removed comment beginnings in Log sections of .cxx files
32Suppressed most violations of coding rules
33
34Revision 1.1.2.1 2000/06/07 14:44:53 gosset
35Addition of files for track reconstruction in C++
36*/
37
38//__________________________________________________________________________
39//
40// Segment for reconstruction in ALICE dimuon spectrometer:
41// two hits for reconstruction in the two chambers of one station
42//__________________________________________________________________________
43
44#include "AliMUONSegment.h"
45
46#include "AliMUON.h"
47#include "AliMUONHitForRec.h"
48#include "AliMUONTrackParam.h"
49#include "AliMUONChamber.h"
50#include "AliRun.h"
51
52ClassImp(AliMUONSegment) // Class implementation in ROOT context
53
54 //__________________________________________________________________________
55AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
56{
57 // Constructor for AliMUONSegment from two HitForRec's,
58 // one, in the first chamber of the station, pointed to by "Hit1",
59 // the other one, in the second chamber of the station, pointed to by "Hit1".
60 // Fills the pointers to both hits,
61 // the slope, the covariance for (coordinate in first chamber, slope),
62 // and the impact parameter at vertex (Z=0),
63 // in bending and non bending planes.
64 // Puts the "fInTrack" flag to "kFALSE".
65 Double_t dz;
66 // pointers to HitForRec's
67 fHitForRecPtr1 = Hit1;
68 fHitForRecPtr2 = Hit2;
69 dz = Hit1->GetZ() - Hit2->GetZ();
70 // bending plane
71 fBendingCoor = Hit1->GetBendingCoor();
72 fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
73 fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
74 fBendingCoorReso2 = Hit1->GetBendingReso2();
75 fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
76 Hit2->GetBendingReso2() ) / dz / dz;
77 fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
78 // non bending plane
79 fNonBendingCoor = Hit1->GetNonBendingCoor();
80 fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
81 fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
82 fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
83 fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
84 Hit2->GetNonBendingReso2() ) / dz / dz;
85 fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
86 // "fInTrack" flag to "kFALSE"
87 fInTrack = kFALSE;
88 return;
89}
90
91AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
92{
93// Dummy copy constructor
94}
95
96AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
97{
98// Dummy assignment operator
99 return *this;
100}
101
a9e2aefa 102 //__________________________________________________________________________
103Int_t AliMUONSegment::Compare(TObject* Segment)
104{
105 // "Compare" function to sort with increasing absolute value
106 // of the "impact parameter" in bending plane.
107 // Returns -1 (0, +1) if |impact parameter| of current Segment
108 // is smaller than (equal to, larger than) |impact parameter| of Segment
109 if (TMath::Abs(((AliMUONSegment*)this)->fBendingSlope)
110 < TMath::Abs(((AliMUONSegment*)Segment)->fBendingSlope))
111 return(-1);
112 // continuous parameter, hence no need for testing equal case
113 else return(+1);
114}
115
116 //__________________________________________________________________________
117Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut)
118{
119 // Calculate the normalized Chi2 between the current Segment (this)
120 // and the Segment pointed to by "Segment",
121 // i.e. the square deviations between the coordinates and the slopes,
122 // in both the bending and the non bending plane,
123 // divided by the variance of the same quantities and by "Sigma2Cut".
124 // Returns 5 if none of the 4 quantities is OK,
125 // something smaller than or equal to 4 otherwise.
126 // Would it be more correct to use a real chi square
127 // including the non diagonal term ????
128 Double_t chi2, chi2Max, diff, normDiff;
129 chi2 = 0.0;
130 chi2Max = 5.0;
131 // coordinate in bending plane
132 diff = this->fBendingCoor - Segment->fBendingCoor;
133 normDiff = diff * diff /
134 (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
135 if (normDiff > 1.0) return chi2Max;
136 chi2 = chi2 + normDiff;
137 // slope in bending plane
138 diff = this->fBendingSlope - Segment->fBendingSlope;
139 normDiff = diff * diff /
140 (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
141 if (normDiff > 1.0) return chi2Max;
142 chi2 = chi2 + normDiff;
143 // coordinate in non bending plane
144 diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
145 normDiff = diff * diff /
146 (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
147 if (normDiff > 1.0) return chi2Max;
148 chi2 = chi2 + normDiff;
149 // slope in non bending plane
150 diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
151 normDiff = diff * diff /
152 (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
153 if (normDiff > 1.0) return chi2Max;
154 chi2 = chi2 + normDiff;
155 return chi2;
156}
157
158 //__________________________________________________________________________
159AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor)
160{
161 // Extrapolates linearly the current Segment (this) to station (0..) "Station".
162 // Multiple Coulomb scattering calculated from "MCSfactor"
163 // corresponding to one chamber,
164 // with one chamber for the coordinate, two chambers for the angle,
165 // due to the arrangement in stations.
166 // Valid from station(1..) 4 to 5 or vice versa.
167 // Returns the pointer to the created AliMUONSegment object
168 // corresponding to this extrapolation.
169 // The caller has the responsibility to delete this object.
170 AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
171 // dZ from first hit of current Segment to first chamber of station "Station"
172 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
173 Double_t dZ =
174 (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
175 // Data in bending plane
176 // coordinate
177 extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
178 // slope
179 extrapSegment->fBendingSlope = this->fBendingSlope;
180 // covariance, including multiple Coulomb scattering over dZ due to one chamber
181 extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
182 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
183 extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
184 extrapSegment->fBendingCoorSlopeReso2 =
185 this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
186 // Data in non bending plane
187 // coordinate
188 extrapSegment->fNonBendingCoor =
189 this->fNonBendingCoor + this->fNonBendingSlope * dZ;
190 // slope
191 extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
192 // covariance, including multiple Coulomb scattering over dZ due to one chamber
193 extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
194 (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
195 extrapSegment->fNonBendingSlopeReso2 =
196 this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
197 extrapSegment->fNonBendingCoorSlopeReso2 =
198 this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
199 return extrapSegment;
200}
201
202 //__________________________________________________________________________
203AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor)
204{
205 // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
206 // Multiple Coulomb scattering calculated from "MCSfactor"
207 // corresponding to one chamber.
208 // Valid from station(1..) 4 to 5 or vice versa.
209 // Returns the pointer to the created AliMUONHitForRec object
210 // corresponding to this extrapolation.
211 // The caller has the responsibility to delete this object.
212 AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
213 // dZ from first hit of current Segment to chamber
214 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
215 Double_t dZ =
216 (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
217 // Data in bending plane
218 // coordinate
219 extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
220 // covariance, including multiple Coulomb scattering over dZ due to one chamber
221 extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
222 (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
223 // Data in non bending plane
224 // coordinate
225 extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
226 this->fNonBendingSlope * dZ);
227 // covariance, including multiple Coulomb scattering over dZ due to one chamber
228 extrapHitForRec->
229 SetNonBendingReso2(this->fNonBendingCoorReso2 +
230 (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
231 return extrapHitForRec;
232}
233
234 //__________________________________________________________________________
04b5ea16 235void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
a9e2aefa 236{
237 // Fill data members with values calculated from the array of track parameters
238 // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
239 // of the station, respectively).
240 // Multiple Coulomb scattering is taking into account with "MCSfactor"
241 // corresponding to one chamber,
242 // with one chamber for the coordinate, two chambers for the angle,
243 // due to the arrangement in stations.
244 // Resolution coming from:
04b5ea16 245 // coordinate in closest station at "Dz1" from current "Station",
a9e2aefa 246 // slope between closest stations, with "Dz2" interval between them,
04b5ea16 247 // interval "Dz3" between chambers of closest station,
248 // extrapolation over "Dz1" from closest station,
249 // "InverseMomentum".
a9e2aefa 250 // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
251 // are assumed to be filled
252 // with the variance on bending and non bending coordinates.
04b5ea16 253 // The "road" is parametrized from the old reco_muon.F
254 // with 8 cm between stations.
a9e2aefa 255 AliMUONTrackParam *param0;
256 Double_t cReso2, sReso2;
04b5ea16 257 // parameters to define the widths of the searching roads in station 0,1,2
258 // width = p0 + p1/ (momentum)^2
259 // station number: 0 1 2
260 static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
261 static Double_t p1BendingCoor[3] = { 986., 821., 446. };
262 static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
263 static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
264 static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
265 static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
266 static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
267 static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
a9e2aefa 268 param0 = &(TrackParam[0]);
04b5ea16 269
270// OLD version
271// // Bending plane
272// fBendingCoor = param0->GetBendingCoor(); // coordinate
273// fBendingSlope = param0->GetBendingSlope(); // slope
274// cReso2 = fBendingCoorReso2;
275// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
276// fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
277// fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
278// // Non bending plane
279// fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
280// fNonBendingSlope = param0->GetNonBendingSlope(); // slope
281// cReso2 = fNonBendingCoorReso2;
282// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
283// fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
284// fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
285
286 // Coordinate and slope
a9e2aefa 287 // Bending plane
288 fBendingCoor = param0->GetBendingCoor(); // coordinate
289 fBendingSlope = param0->GetBendingSlope(); // slope
a9e2aefa 290 // Non bending plane
291 fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
292 fNonBendingSlope = param0->GetNonBendingSlope(); // slope
04b5ea16 293
294 // Resolutions
295 // cReso2 and sReso2 have to be subtracted here from the parametrization
296 // because they are added in the functions "NormalizedChi2WithSegment"
297 // and "NormalizedChi2WithHitForRec"
298 // Bending plane
299 cReso2 = fBendingCoorReso2;
300 sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
301 fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
302 fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
303 // Non bending plane
a9e2aefa 304 cReso2 = fNonBendingCoorReso2;
04b5ea16 305 sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
306 fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
307 fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
a9e2aefa 308 return;
309}
04b5ea16 310
311// OLD function, with roads automatically calculated instead from being parametrized
312// kept because it would be a better solution,
313// if one can really find the right values.
314// //__________________________________________________________________________
315// void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
316// {
317// // Fill data members with values calculated from the array of track parameters
318// // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
319// // of the station, respectively).
320// // Multiple Coulomb scattering is taking into account with "MCSfactor"
321// // corresponding to one chamber,
322// // with one chamber for the coordinate, two chambers for the angle,
323// // due to the arrangement in stations.
324// // Resolution coming from:
325// // coordinate in closest station at "Dz1",
326// // slope between closest stations, with "Dz2" interval between them,
327// // extrapolation over "Dz" from closest station.
328// // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
329// // are assumed to be filled
330// // with the variance on bending and non bending coordinates.
331// AliMUONTrackParam *param0;
332// Double_t cReso2, sReso2;
333// param0 = &(TrackParam[0]);
334// // Bending plane
335// fBendingCoor = param0->GetBendingCoor(); // coordinate
336// fBendingSlope = param0->GetBendingSlope(); // slope
337// cReso2 = fBendingCoorReso2;
338// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
339// fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
340// fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
341// // Non bending plane
342// fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
343// fNonBendingSlope = param0->GetNonBendingSlope(); // slope
344// cReso2 = fNonBendingCoorReso2;
345// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
346// fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
347// fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
348// return;
349// }