Code from MUON-dev joined
[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.1.2.4  2000/06/12 10:10:21  morsch
19 Dummy copy constructor and assignment operator added
20
21 Revision 1.1.2.3  2000/06/09 21:01:16  morsch
22 Make includes consistent with new file structure.
23
24 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
25 Removed comment beginnings in Log sections of .cxx files
26 Suppressed most violations of coding rules
27
28 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
29 Addition of files for track reconstruction in C++
30 */
31
32 //__________________________________________________________________________
33 //
34 // Segment for reconstruction in ALICE dimuon spectrometer:
35 // two hits for reconstruction in the two chambers of one station
36 //__________________________________________________________________________
37
38 #include "AliMUONSegment.h" 
39
40 #include "AliMUON.h"
41 #include "AliMUONHitForRec.h" 
42 #include "AliMUONTrackParam.h" 
43 #include "AliMUONChamber.h" 
44 #include "AliRun.h"
45
46 ClassImp(AliMUONSegment) // Class implementation in ROOT context
47
48   //__________________________________________________________________________
49 AliMUONSegment::AliMUONSegment(AliMUONHitForRec* Hit1, AliMUONHitForRec* Hit2)
50 {
51   // Constructor for AliMUONSegment from two HitForRec's,
52   // one, in the first chamber of the station, pointed to by "Hit1",
53   // the other one, in the second chamber of the station, pointed to by "Hit1".
54   // Fills the pointers to both hits,
55   // the slope, the covariance for (coordinate in first chamber, slope),
56   // and the impact parameter at vertex (Z=0),
57   // in bending and non bending planes.
58   // Puts the "fInTrack" flag to "kFALSE".
59   Double_t dz;
60   // pointers to HitForRec's
61   fHitForRecPtr1 = Hit1;
62   fHitForRecPtr2 = Hit2;
63   dz = Hit1->GetZ() - Hit2->GetZ();
64   // bending plane
65   fBendingCoor = Hit1->GetBendingCoor();
66   fBendingSlope = (fBendingCoor - Hit2->GetBendingCoor()) / dz;
67   fBendingImpact = fBendingCoor - Hit1->GetZ() * fBendingSlope;
68   fBendingCoorReso2 = Hit1->GetBendingReso2();
69   fBendingSlopeReso2 = ( Hit1->GetBendingReso2() +
70                          Hit2->GetBendingReso2() ) / dz / dz;
71   fBendingCoorSlopeReso2 = Hit1->GetBendingReso2() / dz;
72   // non bending plane
73   fNonBendingCoor = Hit1->GetNonBendingCoor();
74   fNonBendingSlope = (fNonBendingCoor - Hit2->GetNonBendingCoor()) / dz;
75   fNonBendingImpact = fNonBendingCoor - Hit1->GetZ() * fNonBendingSlope;
76   fNonBendingCoorReso2 = Hit1->GetNonBendingReso2();
77   fNonBendingSlopeReso2 = ( Hit1->GetNonBendingReso2() +
78                             Hit2->GetNonBendingReso2() ) / dz / dz;
79   fNonBendingCoorSlopeReso2 = Hit1->GetNonBendingReso2() / dz;
80   // "fInTrack" flag to "kFALSE"
81   fInTrack = kFALSE;
82   return;
83 }
84
85 AliMUONSegment::AliMUONSegment (const AliMUONSegment& MUONSegment)
86 {
87 // Dummy copy constructor
88 }
89
90 AliMUONSegment & AliMUONSegment::operator=(const AliMUONSegment& MUONSegment)
91 {
92 // Dummy assignment operator
93     return *this;
94 }
95
96 // Inline functions for Get and Set
97 inline AliMUONHitForRec* AliMUONSegment::GetHitForRec1(void) {
98   // Get fHitForRecPtr1
99   return fHitForRecPtr1;}
100 inline AliMUONHitForRec* AliMUONSegment::GetHitForRec2(void) {
101   // Get fHitForRecPtr2
102   return fHitForRecPtr2;}
103 inline Double_t AliMUONSegment::GetBendingCoorReso2(void) {
104   // Get fBendingCoorReso2
105   return fBendingCoorReso2;}
106 inline void AliMUONSegment::SetBendingCoorReso2(Double_t BendingCoorReso2) {
107   // Set fBendingCoorReso2
108   fBendingCoorReso2 = BendingCoorReso2;}
109 inline Double_t AliMUONSegment::GetNonBendingCoorReso2(void) {
110   // Get fNonBendingCoorReso2
111   return fNonBendingCoorReso2;}
112 inline void AliMUONSegment::SetNonBendingCoorReso2(Double_t NonBendingCoorReso2) {
113   // Set fNonBendingCoorReso2
114   fNonBendingCoorReso2 = NonBendingCoorReso2;}
115 inline Double_t AliMUONSegment::GetBendingImpact(void) {
116   // Get fBendingImpact
117   return fBendingImpact;}
118 inline Bool_t AliMUONSegment::GetInTrack(void) {
119   // Get fInTrack
120   return fInTrack;}
121 inline void AliMUONSegment::SetInTrack(Bool_t InTrack) {
122   // Set fInTrack
123   fInTrack = InTrack;}
124
125   //__________________________________________________________________________
126 Int_t AliMUONSegment::Compare(TObject* Segment)
127 {
128   // "Compare" function to sort with increasing absolute value
129   // of the "impact parameter" in bending plane.
130   // Returns -1 (0, +1) if |impact parameter| of current Segment
131   // is smaller than (equal to, larger than) |impact parameter| of Segment
132   if (TMath::Abs(((AliMUONSegment*)this)->fBendingSlope)
133       < TMath::Abs(((AliMUONSegment*)Segment)->fBendingSlope))
134     return(-1);
135   // continuous parameter, hence no need for testing equal case
136   else return(+1);
137 }
138
139   //__________________________________________________________________________
140 Double_t AliMUONSegment::NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut)
141 {
142   // Calculate the normalized Chi2 between the current Segment (this)
143   // and the Segment pointed to by "Segment",
144   // i.e. the square deviations between the coordinates and the slopes,
145   // in both the bending and the non bending plane,
146   // divided by the variance of the same quantities and by "Sigma2Cut".
147   // Returns 5 if none of the 4 quantities is OK,
148   // something smaller than or equal to 4 otherwise.
149   // Would it be more correct to use a real chi square
150   // including the non diagonal term ????
151   Double_t chi2, chi2Max, diff, normDiff;
152   chi2 = 0.0;
153   chi2Max = 5.0;
154   // coordinate in bending plane
155   diff = this->fBendingCoor - Segment->fBendingCoor;
156   normDiff = diff * diff /
157     (this->fBendingCoorReso2 + Segment->fBendingCoorReso2) / Sigma2Cut;
158   if (normDiff > 1.0) return chi2Max;
159   chi2 = chi2 + normDiff;
160   // slope in bending plane
161   diff = this->fBendingSlope - Segment->fBendingSlope;
162   normDiff = diff * diff /
163     (this->fBendingSlopeReso2 + Segment->fBendingSlopeReso2) / Sigma2Cut;
164   if (normDiff > 1.0) return chi2Max;
165   chi2 = chi2 + normDiff;
166   // coordinate in non bending plane
167   diff = this->fNonBendingCoor - Segment->fNonBendingCoor;
168   normDiff = diff * diff /
169     (this->fNonBendingCoorReso2 + Segment->fNonBendingCoorReso2) / Sigma2Cut;
170   if (normDiff > 1.0) return chi2Max;
171   chi2 = chi2 + normDiff;
172   // slope in non bending plane
173   diff = this->fNonBendingSlope - Segment->fNonBendingSlope;
174   normDiff = diff * diff /
175     (this->fNonBendingSlopeReso2 + Segment->fNonBendingSlopeReso2) / Sigma2Cut;
176   if (normDiff > 1.0) return chi2Max;
177   chi2 = chi2 + normDiff;
178   return chi2;
179 }
180
181   //__________________________________________________________________________
182 AliMUONSegment* AliMUONSegment::CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor)
183 {
184   // Extrapolates linearly the current Segment (this) to station (0..) "Station".
185   // Multiple Coulomb scattering calculated from "MCSfactor"
186   // corresponding to one chamber,
187   // with one chamber for the coordinate, two chambers for the angle,
188   // due to the arrangement in stations.
189   // Valid from station(1..) 4 to 5 or vice versa.
190   // Returns the pointer to the created AliMUONSegment object
191   // corresponding to this extrapolation.
192   // The caller has the responsibility to delete this object.
193   AliMUONSegment* extrapSegment = new AliMUONSegment(); // creates empty new segment
194   // dZ from first hit of current Segment to first chamber of station "Station"
195   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
196   Double_t dZ =
197     (&(pMUON->Chamber(2 * Station)))->Z() - (this->fHitForRecPtr1)->GetZ();
198   // Data in bending plane
199   //  coordinate
200   extrapSegment->fBendingCoor = this->fBendingCoor + this->fBendingSlope * dZ;
201   //  slope
202   extrapSegment->fBendingSlope = this->fBendingSlope;
203   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
204   extrapSegment->fBendingCoorReso2 = this->fBendingCoorReso2 +
205     (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ; // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
206   extrapSegment->fBendingSlopeReso2 = this->fBendingSlopeReso2 + 2.0 * MCSfactor;
207   extrapSegment->fBendingCoorSlopeReso2 =
208     this->fBendingCoorSlopeReso2 + this->fBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
209   // Data in non bending plane
210   //  coordinate
211   extrapSegment->fNonBendingCoor =
212     this->fNonBendingCoor + this->fNonBendingSlope * dZ;
213   //  slope
214   extrapSegment->fNonBendingSlope = this->fNonBendingSlope;
215   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
216   extrapSegment->fNonBendingCoorReso2 = this->fNonBendingCoorReso2 +
217     (this->fNonBendingSlopeReso2 + MCSfactor) *dZ * dZ; // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
218   extrapSegment->fNonBendingSlopeReso2 =
219     this->fNonBendingSlopeReso2 + 2.0 * MCSfactor;
220   extrapSegment->fNonBendingCoorSlopeReso2 =
221     this->fNonBendingCoorSlopeReso2 + this->fNonBendingSlopeReso2 * dZ; // missing: contribution from multiple Coulomb scattering !!!!
222   return extrapSegment;
223 }
224
225   //__________________________________________________________________________
226 AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor)
227 {
228   // Extrapolates linearly the current Segment (this) to chamber(0..) "Chamber".
229   // Multiple Coulomb scattering calculated from "MCSfactor"
230   // corresponding to one chamber.
231   // Valid from station(1..) 4 to 5 or vice versa.
232   // Returns the pointer to the created AliMUONHitForRec object
233   // corresponding to this extrapolation.
234   // The caller has the responsibility to delete this object.
235   AliMUONHitForRec* extrapHitForRec = new AliMUONHitForRec(); // creates empty new HitForRec
236   // dZ from first hit of current Segment to chamber
237   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
238   Double_t dZ =
239     (&(pMUON->Chamber(Chamber)))->Z() - (this->fHitForRecPtr1)->GetZ();
240   // Data in bending plane
241   //  coordinate
242   extrapHitForRec->SetBendingCoor(this->fBendingCoor + this->fBendingSlope * dZ);
243   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
244   extrapHitForRec->SetBendingReso2(this->fBendingCoorReso2 +
245                                    (this->fBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fBendingCoorSlopeReso2 * dZ" !!!!
246   // Data in non bending plane
247   //  coordinate
248  extrapHitForRec ->SetNonBendingCoor(this->fNonBendingCoor +
249                                      this->fNonBendingSlope * dZ);
250   //  covariance, including multiple Coulomb scattering over dZ due to one chamber
251   extrapHitForRec->
252     SetNonBendingReso2(this->fNonBendingCoorReso2 +
253                        (this->fNonBendingSlopeReso2 + MCSfactor) * dZ * dZ); // missing non diagonal term: "2.0 * this->fNonBendingCoorSlopeReso2 * dZ" !!!!
254   return extrapHitForRec;
255 }
256
257   //__________________________________________________________________________
258 void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
259 {
260   // Fill data members with values calculated from the array of track parameters
261   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
262   // of the station, respectively).
263   // Multiple Coulomb scattering is taking into account with "MCSfactor"
264   // corresponding to one chamber,
265   // with one chamber for the coordinate, two chambers for the angle,
266   // due to the arrangement in stations.
267   // Resolution coming from:
268   // coordinate in closest station at "Dz1",
269   // slope between closest stations, with "Dz2" interval between them,
270   // extrapolation over "Dz" from closest station.
271   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
272   // are assumed to be filled
273   // with the variance on bending and non bending coordinates.
274   AliMUONTrackParam *param0;
275   Double_t cReso2, sReso2;
276   param0 = &(TrackParam[0]);
277   // Bending plane
278   fBendingCoor = param0->GetBendingCoor(); // coordinate
279   fBendingSlope = param0->GetBendingSlope(); // slope
280   cReso2 = fBendingCoorReso2;
281   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
282   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
283   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
284   // Non bending plane
285   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
286   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
287   cReso2 = fNonBendingCoorReso2;
288   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
289   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
290   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
291   return;
292 }