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