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