01671570e80cc2e02eb63c71086850806c7e275b
[u/mrichter/AliRoot.git] / MUON / AliMUONSegment.cxx
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$
18 Revision 1.3  2000/06/25 13:06:39  hristov
19 Inline functions moved from *.cxx to *.h files instead of forward declarations
20
21 Revision 1.2  2000/06/15 07:58:48  morsch
22 Code from MUON-dev joined
23
24 Revision 1.1.2.4  2000/06/12 10:10:21  morsch
25 Dummy copy constructor and assignment operator added
26
27 Revision 1.1.2.3  2000/06/09 21:01:16  morsch
28 Make includes consistent with new file structure.
29
30 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
31 Removed comment beginnings in Log sections of .cxx files
32 Suppressed most violations of coding rules
33
34 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
35 Addition 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
52 ClassImp(AliMUONSegment) // Class implementation in ROOT context
53
54   //__________________________________________________________________________
55 AliMUONSegment::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
91 AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
92 {
93 // Dummy copy constructor
94 }
95
96 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
97 {
98 // Dummy assignment operator
99     return *this;
100 }
101
102   //__________________________________________________________________________
103 Int_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   //__________________________________________________________________________
117 Double_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   //__________________________________________________________________________
159 AliMUONSegment* 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   //__________________________________________________________________________
203 AliMUONHitForRec* 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   //__________________________________________________________________________
235 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
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:
245   // coordinate in closest station at "Dz1" from current "Station",
246   // slope between closest stations, with "Dz2" interval between them,
247   // interval "Dz3" between chambers of closest station,
248   // extrapolation over "Dz1" from closest station,
249   // "InverseMomentum".
250   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
251   // are assumed to be filled
252   // with the variance on bending and non bending coordinates.
253   // The "road" is parametrized from the old reco_muon.F
254   // with 8 cm between stations.
255   AliMUONTrackParam *param0;
256   Double_t cReso2, sReso2;
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. };  
268   param0 = &(TrackParam[0]);
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
287   // Bending plane
288   fBendingCoor = param0->GetBendingCoor(); // coordinate
289   fBendingSlope = param0->GetBendingSlope(); // slope
290   // Non bending plane
291   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
292   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
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
304   cReso2 = fNonBendingCoorReso2;
305   sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
306   fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
307   fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
308   return;
309 }
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 // }