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$ |
f6e92cf4 |
18 | Revision 1.2 2000/06/15 07:58:48 morsch |
19 | Code from MUON-dev joined |
20 | |
a9e2aefa |
21 | Revision 1.1.2.3 2000/06/12 10:11:34 morsch |
22 | Dummy copy constructor and assignment operator added |
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 | // Reconstructed track in ALICE dimuon spectrometer |
35 | //__________________________________________________________________________ |
36 | |
37 | #include "AliMUONTrack.h" |
38 | |
39 | #include <iostream.h> |
40 | |
41 | #include <TClonesArray.h> |
42 | #include <TMinuit.h> |
43 | |
44 | #include "AliMUONEventReconstructor.h" |
45 | #include "AliMUONHitForRec.h" |
46 | #include "AliMUONSegment.h" |
47 | #include "AliMUONTrackHit.h" |
48 | |
f6e92cf4 |
49 | #include <stdlib.h> |
50 | |
a9e2aefa |
51 | static AliMUONTrack *trackBeingFitted; |
52 | static AliMUONTrackParam *trackParamBeingFitted; |
53 | |
54 | // Function to be minimized with Minuit |
55 | void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); |
56 | |
57 | ClassImp(AliMUONTrack) // Class implementation in ROOT context |
58 | |
59 | //__________________________________________________________________________ |
60 | AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor) |
61 | { |
62 | // Constructor from two Segment's |
63 | fEventReconstructor = EventReconstructor; // link back to EventReconstructor |
64 | // memory allocation for the TClonesArray of reconstructed TrackHit's |
65 | fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10); |
66 | fNTrackHits = 0; |
67 | AddSegment(BegSegment); // add hits from BegSegment |
68 | AddSegment(EndSegment); // add hits from EndSegment |
69 | fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z |
70 | SetTrackParamAtVertex(); // set track parameters at vertex |
71 | return; |
72 | } |
73 | |
74 | //__________________________________________________________________________ |
75 | AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor) |
76 | { |
77 | // Constructor from one Segment and one HitForRec |
78 | fEventReconstructor = EventReconstructor; // link back to EventReconstructor |
79 | // memory allocation for the TClonesArray of reconstructed TrackHit's |
80 | fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10); |
81 | fNTrackHits = 0; |
82 | AddSegment(Segment); // add hits from Segment |
83 | AddHitForRec(HitForRec); // add HitForRec |
84 | fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z |
85 | SetTrackParamAtVertex(); // set track parameters at vertex |
86 | return; |
87 | } |
88 | |
89 | AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack) |
90 | { |
91 | // Dummy copy constructor |
92 | } |
93 | |
94 | AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack) |
95 | { |
96 | // Dummy assignment operator |
97 | return *this; |
98 | } |
99 | |
100 | // Inline functions for Get and Set: inline removed because it does not work !!!! |
101 | AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) { |
102 | // Get pointer to fTrackParamAtVertex |
103 | return &fTrackParamAtVertex;} |
104 | AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) { |
105 | // Get pointer to TrackParamAtFirstHit |
106 | return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} |
107 | TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) { |
108 | // Get fTrackHitsPtr |
109 | return fTrackHitsPtr;} |
110 | Int_t AliMUONTrack::GetNTrackHits(void) { |
111 | // Get fNTrackHits |
112 | return fNTrackHits;} |
113 | |
114 | //__________________________________________________________________________ |
115 | void AliMUONTrack::RecursiveDump(void) |
116 | { |
117 | // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's |
118 | AliMUONTrackHit *trackHit; |
119 | AliMUONHitForRec *hitForRec; |
120 | cout << "Recursive dump of Track: " << this << endl; |
121 | // Track |
122 | this->Dump(); |
123 | for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) { |
124 | trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]); |
125 | // TrackHit |
126 | cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl; |
127 | trackHit->Dump(); |
128 | hitForRec = trackHit->GetHitForRecPtr(); |
129 | // HitForRec |
130 | cout << "HitForRec: " << hitForRec << endl; |
131 | hitForRec->Dump(); |
132 | } |
133 | return; |
134 | } |
135 | |
136 | //__________________________________________________________________________ |
137 | void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam) |
138 | { |
139 | // Fit the current track ("this"), |
140 | // starting with track parameters pointed to by "TrackParam", |
141 | // and with 3 or 5 parameters ("NParam"): |
142 | // 3 if one keeps X and Y fixed in "TrackParam", |
143 | // 5 if one lets them vary. |
144 | if ((NParam != 3) && (NParam != 5)) { |
145 | cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam; |
146 | cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl; |
147 | exit(0); // right instruction for exit ???? |
148 | } |
149 | Int_t error = 0; |
150 | Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y; |
151 | TString parName; |
152 | TMinuit *minuit = new TMinuit(5); |
153 | trackBeingFitted = this; // for the track to be known from the function to minimize |
154 | trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ???? |
155 | minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ???? |
156 | // how to go faster ???? choice of Minuit parameters like EDM ???? |
157 | minuit->SetFCN(TrackChi2); |
158 | minuit->SetPrintLevel(1); // More printing !!!! |
159 | // set first 3 parameters (try with no limits for the search: min=max=0) |
160 | minuit->mnparm(0, "InvBenP", |
161 | TrackParam->GetInverseBendingMomentum(), |
162 | 0.003, 0.0, 0.0, error); |
163 | minuit->mnparm(1, "BenS", |
164 | TrackParam->GetBendingSlope(), |
165 | 0.001, 0.0, 0.0, error); |
166 | minuit->mnparm(2, "NonBenS", |
167 | TrackParam->GetNonBendingSlope(), |
168 | 0.001, 0.0, 0.0, error); |
169 | if (NParam == 5) { |
170 | // set last 2 parameters (try with no limits for the search: min=max=0) |
171 | minuit->mnparm(3, "X", |
172 | TrackParam->GetNonBendingCoor(), |
173 | 0.03, 0.0, 0.0, error); |
174 | minuit->mnparm(4, "Y", |
175 | TrackParam->GetBendingCoor(), |
176 | 0.10, 0.0, 0.0, error); |
177 | } |
178 | // search without gradient calculation in the function |
179 | minuit->mnexcm("SET NOGRADIENT", arg, 0, error); |
180 | // minimization |
181 | minuit->mnexcm("MINIMIZE", arg, 0, error); |
182 | // exit from Minuit |
183 | minuit->mnexcm("EXIT", arg, 0, error); // necessary ???? |
184 | // print results |
185 | minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error); |
186 | minuit->mnpout(1, parName, benC, errorParam, lower, upper, error); |
187 | minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error); |
188 | if (NParam == 5) { |
189 | minuit->mnpout(3, parName, x, errorParam, lower, upper, error); |
190 | minuit->mnpout(4, parName, y, errorParam, lower, upper, error); |
191 | } |
192 | // result of the fit into track parameters |
193 | TrackParam->SetInverseBendingMomentum(invBenP); |
194 | TrackParam->SetBendingSlope(benC); |
195 | TrackParam->SetNonBendingSlope(nonBenC); |
196 | if (NParam == 5) { |
197 | TrackParam->SetNonBendingCoor(x); |
198 | TrackParam->SetBendingCoor(y); |
199 | } |
200 | trackBeingFitted = NULL; |
201 | delete minuit; |
202 | } |
203 | |
204 | //__________________________________________________________________________ |
205 | void AliMUONTrack::AddSegment(AliMUONSegment* Segment) |
206 | { |
207 | // Add Segment |
208 | AddHitForRec(Segment->GetHitForRec1()); // 1st hit |
209 | AddHitForRec(Segment->GetHitForRec2()); // 2nd hit |
210 | } |
211 | |
212 | //__________________________________________________________________________ |
213 | void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) |
214 | { |
215 | // Add HitForRec |
216 | new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec); |
217 | fNTrackHits++; |
218 | } |
219 | |
220 | //__________________________________________________________________________ |
221 | void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) |
222 | { |
223 | // Set track parameters at TrackHit with index "indexHit" |
224 | // from the track parameters pointed to by "TrackParam". |
225 | AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); |
226 | trackHit->SetTrackParam(TrackParam); |
227 | } |
228 | |
229 | //__________________________________________________________________________ |
230 | void AliMUONTrack::SetTrackParamAtVertex() |
231 | { |
232 | // Set track parameters at vertex. |
233 | // TrackHit's are assumed to be only in stations(1..) 4 and 5, |
234 | // and sorted according to increasing Z.. |
235 | // Parameters are calculated from information in HitForRec's |
236 | // of first and last TrackHit's. |
237 | AliMUONTrackParam *trackParam = |
238 | &fTrackParamAtVertex; // pointer to track parameters |
239 | // Pointer to HitForRec of first TrackHit |
240 | AliMUONHitForRec *firstHit = |
241 | ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr(); |
242 | // Pointer to HitForRec of last TrackHit |
243 | AliMUONHitForRec *lastHit = |
244 | ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr(); |
245 | // Z difference between first and last hits |
246 | Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ(); |
247 | // bending slope in stations(1..) 4 and 5 |
248 | Double_t bendingSlope = |
249 | (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ; |
250 | trackParam->SetBendingSlope(bendingSlope); |
251 | // impact parameter |
252 | Double_t impactParam = |
253 | firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and lastHit ???? |
254 | // signed bending momentum |
255 | Double_t signedBendingMomentum = |
256 | fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam); |
257 | trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum); |
258 | // bending slope at vertex |
259 | trackParam-> |
260 | SetBendingSlope(bendingSlope + |
261 | impactParam / fEventReconstructor->GetSimpleBPosition()); |
262 | // non bending slope |
263 | Double_t nonBendingSlope = |
264 | (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ; |
265 | trackParam->SetNonBendingSlope(nonBendingSlope); |
266 | // vertex coordinates at (0,0,0) |
267 | trackParam->SetZ(0.0); |
268 | trackParam->SetBendingCoor(0.0); |
269 | trackParam->SetNonBendingCoor(0.0); |
270 | } |
271 | |
272 | //__________________________________________________________________________ |
273 | void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag) |
274 | { |
275 | // Return the "Chi2" to be minimized with Minuit for track fitting, |
276 | // with "NParam" parameters |
277 | // and their current values in array pointed to by "Param". |
278 | // Assumes that the track hits are sorted according to increasing Z. |
279 | // Track parameters at each TrackHit are updated accordingly. |
280 | AliMUONTrackHit* hit; |
281 | AliMUONTrackParam param1; |
282 | Int_t hitNumber; |
283 | Double_t zHit; |
284 | Chi2 = 0.0; // initialize Chi2 |
285 | // copy of track parameters to be fitted |
286 | param1 = *trackParamBeingFitted; |
287 | // Minuit parameters to be fitted into this copy |
288 | param1.SetInverseBendingMomentum(Param[0]); |
289 | param1.SetBendingSlope(Param[1]); |
290 | param1.SetNonBendingSlope(Param[2]); |
291 | if (NParam == 5) { |
292 | param1.SetNonBendingCoor(Param[3]); |
293 | param1.SetBendingCoor(Param[4]); |
294 | } |
295 | // Follow track through all planes of track hits |
296 | for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) { |
297 | hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; |
298 | zHit = hit->GetHitForRecPtr()->GetZ(); |
299 | // do something special if 2 hits with same Z ???? |
300 | // security against infinite loop ???? |
301 | (¶m1)->ExtrapToZ(zHit); // extrapolation |
302 | hit->SetTrackParam(¶m1); |
303 | // Increment Chi2 |
304 | // done hit per hit, with hit resolution, |
305 | // and not with point and angle like in "reco_muon.F" !!!! |
306 | // Needs to add multiple scattering contribution ???? |
307 | Double_t dX = |
308 | hit->GetHitForRecPtr()->GetNonBendingCoor() - (¶m1)->GetNonBendingCoor(); |
309 | Double_t dY = |
310 | hit->GetHitForRecPtr()->GetBendingCoor() - (¶m1)->GetBendingCoor(); |
311 | Chi2 = |
312 | Chi2 + |
313 | dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() + |
314 | dY * dY / hit->GetHitForRecPtr()->GetBendingReso2(); |
315 | } |
316 | } |