Modifications used for addendum to Dimuon TDR (JP Cussonneau):
[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.4  2000/06/30 10:15:48  gosset
19 Changes to EventReconstructor...:
20 precision fit with multiple Coulomb scattering;
21 extrapolation to vertex with Branson correction in absorber (JPC)
22
23 Revision 1.3  2000/06/25 13:06:39  hristov
24 Inline functions moved from *.cxx to *.h files instead of forward declarations
25
26 Revision 1.2  2000/06/15 07:58:48  morsch
27 Code from MUON-dev joined
28
29 Revision 1.1.2.4  2000/06/12 10:10:21  morsch
30 Dummy copy constructor and assignment operator added
31
32 Revision 1.1.2.3  2000/06/09 21:01:16  morsch
33 Make includes consistent with new file structure.
34
35 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
36 Removed comment beginnings in Log sections of .cxx files
37 Suppressed most violations of coding rules
38
39 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
40 Addition 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
57 ClassImp(AliMUONSegment) // Class implementation in ROOT context
58
59   //__________________________________________________________________________
60 AliMUONSegment::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
96 AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
97 {
98 // Dummy copy constructor
99 }
100
101 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
102 {
103 // Dummy assignment operator
104     return *this;
105 }
106
107   //__________________________________________________________________________
108 Int_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   //__________________________________________________________________________
122 Double_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   //__________________________________________________________________________
164 AliMUONSegment* 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   //__________________________________________________________________________
208 AliMUONHitForRec* 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   //__________________________________________________________________________
240 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
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:
250   // coordinate in closest station at "Dz1" from current "Station",
251   // slope between closest stations, with "Dz2" interval between them,
252   // interval "Dz3" between chambers of closest station,
253   // extrapolation over "Dz1" from closest station,
254   // "InverseMomentum".
255   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
256   // are assumed to be filled
257   // with the variance on bending and non bending coordinates.
258   // The "road" is parametrized from the old reco_muon.F
259   // with 8 cm between stations.
260   AliMUONTrackParam *param0;
261 //   Double_t cReso2, sReso2;
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
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. };  
282   param0 = &(TrackParam[0]);
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
301   // Bending plane
302   fBendingCoor = param0->GetBendingCoor(); // coordinate
303   fBendingSlope = param0->GetBendingSlope(); // slope
304   // Non bending plane
305   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
306   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
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
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;
317   // Non bending plane
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;
322   return;
323 }
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 // }