Macros named
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
CommitLineData
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$
d0bfce8d 18Revision 1.7 2000/09/19 15:50:46 gosset
19TrackChi2MCS function: covariance matrix better calculated,
20taking into account missing planes...
21
eb9c9dab 22Revision 1.6 2000/07/20 12:45:27 gosset
23New "EventReconstructor..." structure,
24 hopefully more adapted to tree/streamer.
25"AliMUONEventReconstructor::RemoveDoubleTracks"
26 to keep only one track among similar ones.
27
8429a5e4 28Revision 1.5 2000/07/18 16:04:06 gosset
29AliMUONEventReconstructor package:
30* a few minor modifications and more comments
31* a few corrections
32 * right sign for Z of raw clusters
33 * right loop over chambers inside station
34 * symmetrized covariance matrix for measurements (TrackChi2MCS)
35 * right sign of charge in extrapolation (ExtrapToZ)
36 * right zEndAbsorber for Branson correction below 3 degrees
37* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
38* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
39
956019b6 40Revision 1.4 2000/06/30 10:15:48 gosset
41Changes to EventReconstructor...:
42precision fit with multiple Coulomb scattering;
43extrapolation to vertex with Branson correction in absorber (JPC)
44
04b5ea16 45Revision 1.3 2000/06/25 13:23:28 hristov
46stdlib.h needed for non-Linux compilation
47
f6e92cf4 48Revision 1.2 2000/06/15 07:58:48 morsch
49Code from MUON-dev joined
50
a9e2aefa 51Revision 1.1.2.3 2000/06/12 10:11:34 morsch
52Dummy copy constructor and assignment operator added
53
54Revision 1.1.2.2 2000/06/09 12:58:05 gosset
55Removed comment beginnings in Log sections of .cxx files
56Suppressed most violations of coding rules
57
58Revision 1.1.2.1 2000/06/07 14:44:53 gosset
59Addition of files for track reconstruction in C++
60*/
61
62//__________________________________________________________________________
63//
64// Reconstructed track in ALICE dimuon spectrometer
65//__________________________________________________________________________
66
67#include "AliMUONTrack.h"
68
69#include <iostream.h>
70
71#include <TClonesArray.h>
04b5ea16 72#include <TMath.h>
d0bfce8d 73#include <TMatrixD.h>
8429a5e4 74#include <TObjArray.h>
956019b6 75#include <TVirtualFitter.h>
a9e2aefa 76
77#include "AliMUONEventReconstructor.h"
78#include "AliMUONHitForRec.h"
79#include "AliMUONSegment.h"
80#include "AliMUONTrackHit.h"
81
f6e92cf4 82#include <stdlib.h>
83
04b5ea16 84// Functions to be minimized with Minuit
a9e2aefa 85void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
04b5ea16 86void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
87
88Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
a9e2aefa 89
90ClassImp(AliMUONTrack) // Class implementation in ROOT context
91
956019b6 92TVirtualFitter* AliMUONTrack::fgFitter = NULL;
93
a9e2aefa 94 //__________________________________________________________________________
95AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
96{
97 // Constructor from two Segment's
98 fEventReconstructor = EventReconstructor; // link back to EventReconstructor
8429a5e4 99 // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
100 fTrackHitsPtr = new TObjArray(10);
a9e2aefa 101 fNTrackHits = 0;
102 AddSegment(BegSegment); // add hits from BegSegment
103 AddSegment(EndSegment); // add hits from EndSegment
104 fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
105 SetTrackParamAtVertex(); // set track parameters at vertex
8429a5e4 106 // set fit conditions...
04b5ea16 107 fFitMCS = 0;
956019b6 108 fFitNParam = 3;
109 fFitStart = 1;
8429a5e4 110 fFitFMin = -1.0;
a9e2aefa 111 return;
112}
113
114 //__________________________________________________________________________
115AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor)
116{
117 // Constructor from one Segment and one HitForRec
118 fEventReconstructor = EventReconstructor; // link back to EventReconstructor
8429a5e4 119 // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
120 fTrackHitsPtr = new TObjArray(10);
a9e2aefa 121 fNTrackHits = 0;
122 AddSegment(Segment); // add hits from Segment
123 AddHitForRec(HitForRec); // add HitForRec
124 fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
125 SetTrackParamAtVertex(); // set track parameters at vertex
8429a5e4 126 // set fit conditions...
04b5ea16 127 fFitMCS = 0;
956019b6 128 fFitNParam = 3;
129 fFitStart = 1;
8429a5e4 130 fFitFMin = -1.0;
a9e2aefa 131 return;
132}
133
8429a5e4 134 //__________________________________________________________________________
135AliMUONTrack::~AliMUONTrack()
136{
137 // Destructor
138 if (fTrackHitsPtr) {
139 delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's
140 fTrackHitsPtr = NULL;
141 }
142}
143
956019b6 144 //__________________________________________________________________________
a9e2aefa 145AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
146{
147// Dummy copy constructor
148}
149
956019b6 150 //__________________________________________________________________________
a9e2aefa 151AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
152{
153// Dummy assignment operator
154 return *this;
155}
156
8429a5e4 157 //__________________________________________________________________________
158void AliMUONTrack::Remove()
159{
160 // Remove current track from array of tracks,
161 // and corresponding track hits from array of track hits.
162 // Compress the TClonesArray it belongs to.
163 AliMUONTrackHit *nextTrackHit;
164 AliMUONEventReconstructor *eventRec = this->fEventReconstructor;
165 TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr();
166 // Loop over all track hits of track
167 AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First();
168 while (trackHit) {
169 nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit);
170 // Remove TrackHit from event TClonesArray.
171 // Destructor is called,
172 // hence links between HitForRec's and TrackHit's are updated
173 trackHitsPtr->Remove(trackHit);
174 trackHit = nextTrackHit;
175 }
176 // Remove the track from event TClonesArray
177 // Destructor is called,
178 // hence space for TObjArray of pointers to TrackHit's is freed
179 eventRec->GetRecTracksPtr()->Remove(this);
180 // Number of tracks decreased by 1
181 eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1);
182 // Compress event TClonesArray of Track's:
183 // this is essential to retrieve the TClonesArray afterwards
184 eventRec->GetRecTracksPtr()->Compress();
185 // Compress event TClonesArray of TrackHit's:
186 // this is probably also essential to retrieve the TClonesArray afterwards
187 trackHitsPtr->Compress();
188}
189
04b5ea16 190 //__________________________________________________________________________
191void AliMUONTrack::SetFitMCS(Int_t FitMCS)
192{
956019b6 193 // Set multiple Coulomb scattering option for track fit "fFitMCS"
04b5ea16 194 // from "FitMCS" argument: 0 without, 1 with
956019b6 195 if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
04b5ea16 196 // better implementation with enum(with, without) ????
956019b6 197 else {
198 cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
199 cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
200 exit(0);
201 }
202 return;
203}
204
205 //__________________________________________________________________________
206void AliMUONTrack::SetFitNParam(Int_t FitNParam)
207{
208 // Set number of parameters for track fit "fFitNParam" from "FitNParam":
209 // 3 for momentum, 5 for momentum and position
210 if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
211 else {
212 cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
213 cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
214 exit(0);
215 }
216 return;
217}
218
219 //__________________________________________________________________________
220void AliMUONTrack::SetFitStart(Int_t FitStart)
221{
222 // Set multiple Coulomb scattering option for track fit "fFitStart"
223 // from "FitStart" argument: 0 without, 1 with
224 if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
225 // better implementation with enum(vertex, firstHit) ????
226 else {
227 cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
228 cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
229 exit(0);
230 }
04b5ea16 231 return;
232}
233
956019b6 234 //__________________________________________________________________________
a9e2aefa 235AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
236 // Get pointer to TrackParamAtFirstHit
237 return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
a9e2aefa 238
239 //__________________________________________________________________________
240void AliMUONTrack::RecursiveDump(void)
241{
242 // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
243 AliMUONTrackHit *trackHit;
244 AliMUONHitForRec *hitForRec;
245 cout << "Recursive dump of Track: " << this << endl;
246 // Track
247 this->Dump();
248 for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
249 trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]);
250 // TrackHit
251 cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
252 trackHit->Dump();
253 hitForRec = trackHit->GetHitForRecPtr();
254 // HitForRec
255 cout << "HitForRec: " << hitForRec << endl;
256 hitForRec->Dump();
257 }
258 return;
259}
260
8429a5e4 261 //__________________________________________________________________________
262Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track)
263{
264 // Returns the number of hits in common
265 // between the current track ("this")
266 // and the track pointed to by "Track".
267 Int_t hitsInCommon = 0;
268 AliMUONTrackHit *trackHit1, *trackHit2;
269 // Loop over hits of first track
270 trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First();
271 while (trackHit1) {
272 // Loop over hits of second track
273 trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First();
274 while (trackHit2) {
275 // Increment "hitsInCommon" if both TrackHits point to the same HitForRec
276 if ( (trackHit1->GetHitForRecPtr()) ==
277 (trackHit2->GetHitForRecPtr()) ) hitsInCommon++;
278 trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2);
279 } // trackHit2
280 trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1);
281 } // trackHit1
282 return hitsInCommon;
283}
284
a9e2aefa 285 //__________________________________________________________________________
956019b6 286void AliMUONTrack::Fit()
a9e2aefa 287{
288 // Fit the current track ("this"),
956019b6 289 // with or without multiple Coulomb scattering according to "fFitMCS",
290 // with the number of parameters given by "fFitNParam"
291 // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
292 // starting, according to "fFitStart",
293 // with track parameters at vertex or at the first TrackHit.
294 // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
295 // by calling the corresponding Set methods.
a9e2aefa 296 Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
956019b6 297 char parName[50];
298 AliMUONTrackParam *trackParam;
299 // Check if Minuit is initialized...
300 fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
301 fgFitter->Clear(); // necessary ???? probably yes
302 // how to reset the printout number at every fit ????
303 // is there any risk to leave it like that ????
a9e2aefa 304 // how to go faster ???? choice of Minuit parameters like EDM ????
04b5ea16 305 // choice of function to be minimized according to fFitMCS
956019b6 306 if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
307 else fgFitter->SetFCN(TrackChi2MCS);
d0bfce8d 308 arg[0] = -1;
956019b6 309 fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
310 // Parameters according to "fFitStart"
311 // (should be a function to be used at every place where needed ????)
312 if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
313 else trackParam = this->GetTrackParamAtFirstHit();
314 // set first 3 Minuit parameters
04b5ea16 315 // could be tried with no limits for the search (min=max=0) ????
956019b6 316 fgFitter->SetParameter(0, "InvBenP",
317 trackParam->GetInverseBendingMomentum(),
318 0.003, -0.4, 0.4);
319 fgFitter->SetParameter(1, "BenS",
320 trackParam->GetBendingSlope(),
321 0.001, -0.5, 0.5);
322 fgFitter->SetParameter(2, "NonBenS",
323 trackParam->GetNonBendingSlope(),
324 0.001, -0.5, 0.5);
325 if (fFitNParam == 5) {
d0bfce8d 326 // set last 2 Minuit parameters
327 // mandatory limits in Bending to avoid NaN values of parameters
956019b6 328 fgFitter->SetParameter(3, "X",
329 trackParam->GetNonBendingCoor(),
d0bfce8d 330 0.03, -500.0, 500.0);
331 // mandatory limits in non Bending to avoid NaN values of parameters
956019b6 332 fgFitter->SetParameter(4, "Y",
333 trackParam->GetBendingCoor(),
d0bfce8d 334 0.10, -500.0, 500.0);
a9e2aefa 335 }
336 // search without gradient calculation in the function
956019b6 337 fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
a9e2aefa 338 // minimization
956019b6 339 fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
a9e2aefa 340 // exit from Minuit
956019b6 341 fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
342 // get results into "invBenP", "benC", "nonBenC" ("x", "y")
343 fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
344 fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
345 fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
346 if (fFitNParam == 5) {
347 fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
348 fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
a9e2aefa 349 }
350 // result of the fit into track parameters
956019b6 351 trackParam->SetInverseBendingMomentum(invBenP);
352 trackParam->SetBendingSlope(benC);
353 trackParam->SetNonBendingSlope(nonBenC);
354 if (fFitNParam == 5) {
355 trackParam->SetNonBendingCoor(x);
356 trackParam->SetBendingCoor(y);
a9e2aefa 357 }
8429a5e4 358 // global result of the fit
359 Double_t fedm, errdef;
360 Int_t npari, nparx;
361 fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx);
a9e2aefa 362}
363
364 //__________________________________________________________________________
365void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
366{
8429a5e4 367 // Add Segment to the track
a9e2aefa 368 AddHitForRec(Segment->GetHitForRec1()); // 1st hit
369 AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
370}
371
372 //__________________________________________________________________________
373void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
374{
8429a5e4 375 // Add HitForRec to the track:
376 // actual TrackHit into TClonesArray of TrackHit's for the event;
377 // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track
378 TClonesArray *recTrackHitsPtr = this->fEventReconstructor->GetRecTrackHitsPtr();
379 Int_t eventTrackHits = this->fEventReconstructor->GetNRecTrackHits();
380 // event
381 AliMUONTrackHit* trackHit =
382 new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec);
383 this->fEventReconstructor->SetNRecTrackHits(eventTrackHits + 1);
384 // track
385 fTrackHitsPtr->Add(trackHit);
a9e2aefa 386 fNTrackHits++;
387}
388
389 //__________________________________________________________________________
390void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
391{
392 // Set track parameters at TrackHit with index "indexHit"
393 // from the track parameters pointed to by "TrackParam".
394 AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
395 trackHit->SetTrackParam(TrackParam);
396}
397
398 //__________________________________________________________________________
399void AliMUONTrack::SetTrackParamAtVertex()
400{
401 // Set track parameters at vertex.
402 // TrackHit's are assumed to be only in stations(1..) 4 and 5,
403 // and sorted according to increasing Z..
404 // Parameters are calculated from information in HitForRec's
405 // of first and last TrackHit's.
406 AliMUONTrackParam *trackParam =
407 &fTrackParamAtVertex; // pointer to track parameters
408 // Pointer to HitForRec of first TrackHit
409 AliMUONHitForRec *firstHit =
410 ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
411 // Pointer to HitForRec of last TrackHit
412 AliMUONHitForRec *lastHit =
413 ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
414 // Z difference between first and last hits
415 Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
416 // bending slope in stations(1..) 4 and 5
417 Double_t bendingSlope =
418 (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
419 trackParam->SetBendingSlope(bendingSlope);
420 // impact parameter
421 Double_t impactParam =
422 firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and lastHit ????
423 // signed bending momentum
424 Double_t signedBendingMomentum =
425 fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam);
426 trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
427 // bending slope at vertex
428 trackParam->
429 SetBendingSlope(bendingSlope +
430 impactParam / fEventReconstructor->GetSimpleBPosition());
431 // non bending slope
432 Double_t nonBendingSlope =
433 (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
434 trackParam->SetNonBendingSlope(nonBendingSlope);
435 // vertex coordinates at (0,0,0)
436 trackParam->SetZ(0.0);
437 trackParam->SetBendingCoor(0.0);
438 trackParam->SetNonBendingCoor(0.0);
439}
440
441 //__________________________________________________________________________
442void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
443{
444 // Return the "Chi2" to be minimized with Minuit for track fitting,
445 // with "NParam" parameters
446 // and their current values in array pointed to by "Param".
447 // Assumes that the track hits are sorted according to increasing Z.
448 // Track parameters at each TrackHit are updated accordingly.
04b5ea16 449 // Multiple Coulomb scattering is not taken into account
956019b6 450 AliMUONTrack *trackBeingFitted;
a9e2aefa 451 AliMUONTrackHit* hit;
452 AliMUONTrackParam param1;
453 Int_t hitNumber;
454 Double_t zHit;
455 Chi2 = 0.0; // initialize Chi2
456 // copy of track parameters to be fitted
956019b6 457 trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
458 if (trackBeingFitted->GetFitStart() == 0)
459 param1 = *(trackBeingFitted->GetTrackParamAtVertex());
460 else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
a9e2aefa 461 // Minuit parameters to be fitted into this copy
462 param1.SetInverseBendingMomentum(Param[0]);
463 param1.SetBendingSlope(Param[1]);
464 param1.SetNonBendingSlope(Param[2]);
465 if (NParam == 5) {
466 param1.SetNonBendingCoor(Param[3]);
467 param1.SetBendingCoor(Param[4]);
468 }
469 // Follow track through all planes of track hits
470 for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
471 hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
472 zHit = hit->GetHitForRecPtr()->GetZ();
473 // do something special if 2 hits with same Z ????
474 // security against infinite loop ????
475 (&param1)->ExtrapToZ(zHit); // extrapolation
476 hit->SetTrackParam(&param1);
477 // Increment Chi2
478 // done hit per hit, with hit resolution,
479 // and not with point and angle like in "reco_muon.F" !!!!
480 // Needs to add multiple scattering contribution ????
481 Double_t dX =
482 hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
483 Double_t dY =
484 hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
485 Chi2 =
486 Chi2 +
487 dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
488 dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
489 }
490}
04b5ea16 491
492 //__________________________________________________________________________
493void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
494{
495 // Return the "Chi2" to be minimized with Minuit for track fitting,
496 // with "NParam" parameters
497 // and their current values in array pointed to by "Param".
498 // Assumes that the track hits are sorted according to increasing Z.
499 // Track parameters at each TrackHit are updated accordingly.
500 // Multiple Coulomb scattering is taken into account with covariance matrix.
956019b6 501 AliMUONTrack *trackBeingFitted;
04b5ea16 502 AliMUONTrackParam param1;
503 Chi2 = 0.0; // initialize Chi2
504 // copy of track parameters to be fitted
956019b6 505 trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
506 if (trackBeingFitted->GetFitStart() == 0)
507 param1 = *(trackBeingFitted->GetTrackParamAtVertex());
508 else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
04b5ea16 509 // Minuit parameters to be fitted into this copy
510 param1.SetInverseBendingMomentum(Param[0]);
511 param1.SetBendingSlope(Param[1]);
512 param1.SetNonBendingSlope(Param[2]);
513 if (NParam == 5) {
514 param1.SetNonBendingCoor(Param[3]);
515 param1.SetBendingCoor(Param[4]);
516 }
517
956019b6 518 AliMUONTrackHit *hit;
519 Bool_t goodDeterminant;
eb9c9dab 520 Int_t chCurrent, chPrev, hitNumber, hitNumber1, hitNumber2, hitNumber3;
521 Double_t z, z1, z2, z3;
522 AliMUONTrackHit *hit1, *hit2, *hit3;
523 Double_t hbc1, hbc2, pbc1, pbc2;
524 Double_t hnbc1, hnbc2, pnbc1, pnbc2;
525 Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
d0bfce8d 526 TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
527 TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
eb9c9dab 528 Double_t *msa2 = new Double_t[numberOfHit];
04b5ea16 529
530 // Predicted coordinates and multiple scattering angles are first calculated
531 for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
532 hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
eb9c9dab 533 z = hit->GetHitForRecPtr()->GetZ();
04b5ea16 534 // do something special if 2 hits with same Z ????
535 // security against infinite loop ????
eb9c9dab 536 (&param1)->ExtrapToZ(z); // extrapolation
04b5ea16 537 hit->SetTrackParam(&param1);
eb9c9dab 538 // square of multiple scattering angle at current hit, with one chamber
539 msa2[hitNumber] = MultipleScatteringAngle2(hit);
540 // correction for eventual missing hits or multiple hits in a chamber,
541 // according to the number of chambers
542 // between the current hit and the previous one
543 chCurrent = hit->GetHitForRecPtr()->GetChamberNumber();
544 if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
545 chPrev = chCurrent;
04b5ea16 546 }
547
548 // Calculates the covariance matrix
eb9c9dab 549 for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
550 hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
551 z1 = hit1->GetHitForRecPtr()->GetZ();
04b5ea16 552 for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
eb9c9dab 553 hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
554 z2 = hit2->GetHitForRecPtr()->GetZ();
956019b6 555 // initialization to 0 (diagonal plus upper triangular part)
556 (*covBending)(hitNumber2, hitNumber1) = 0.0;
557 // contribution from multiple scattering in bending plane:
558 // loop over upstream hits
559 for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {
eb9c9dab 560 hit3 = (AliMUONTrackHit*)
561 (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3];
562 z3 = hit3->GetHitForRecPtr()->GetZ();
04b5ea16 563 (*covBending)(hitNumber2, hitNumber1) =
564 (*covBending)(hitNumber2, hitNumber1) +
eb9c9dab 565 ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]);
956019b6 566 }
567 // equal contribution from multiple scattering in non bending plane
04b5ea16 568 (*covNonBending)(hitNumber2, hitNumber1) =
569 (*covBending)(hitNumber2, hitNumber1);
956019b6 570 if (hitNumber1 == hitNumber2) {
571 // Diagonal elements: add contribution from position measurements
572 // in bending plane
573 (*covBending)(hitNumber2, hitNumber1) =
eb9c9dab 574 (*covBending)(hitNumber2, hitNumber1) +
575 hit1->GetHitForRecPtr()->GetBendingReso2();
956019b6 576 // and in non bending plane
04b5ea16 577 (*covNonBending)(hitNumber2, hitNumber1) =
eb9c9dab 578 (*covNonBending)(hitNumber2, hitNumber1) +
579 hit1->GetHitForRecPtr()->GetNonBendingReso2();
956019b6 580 }
581 else {
582 // Non diagonal elements: symmetrization
583 // for bending plane
584 (*covBending)(hitNumber1, hitNumber2) =
585 (*covBending)(hitNumber2, hitNumber1);
586 // and non bending plane
587 (*covNonBending)(hitNumber1, hitNumber2) =
588 (*covNonBending)(hitNumber2, hitNumber1);
589 }
590 } // for (hitNumber2 = hitNumber1;...
591 } // for (hitNumber1 = 0;...
d0bfce8d 592 // Normalization of covariance matrices
593 Double_t normCovBending2 = covBending->E2Norm();
594 Double_t normCovNonBending2 = covNonBending->E2Norm();
595 (*covBending) *= 1/normCovBending2;
596 (*covNonBending) *= 1/normCovNonBending2;
597// if (covBending->Determinant() < 1.e-33) {
598// printf(" *** covBending *** \n");
599// covBending->Print();
600// printf(" *** covNonBending *** \n");
601// covNonBending->Print();
602// cout << " number of hits " << numberOfHit << endl;
603// cout << "Momentum = " << 1/Param[0] <<endl;
604// cout << "normCovBending = " << normCovBending2 << endl;
605// cout << "normCovNonBending = " << normCovNonBending2 << endl;
606// exit(0);
607
608// }
04b5ea16 609 // Inverts covariance matrix
956019b6 610 goodDeterminant = kTRUE;
611 // check whether the Invert method returns flag if matrix cannot be inverted,
612 // and do not calculate the Determinant in that case !!!!
04b5ea16 613 if (covBending->Determinant() != 0) {
614 covBending->Invert();
615 } else {
956019b6 616 goodDeterminant = kFALSE;
04b5ea16 617 cout << "Warning in ChiMCS Determinant Bending=0: " << endl;
618 }
956019b6 619 if (covNonBending->Determinant() != 0) {
04b5ea16 620 covNonBending->Invert();
621 } else {
956019b6 622 goodDeterminant = kFALSE;
04b5ea16 623 cout << "Warning in ChiMCS Determinant non Bending=0: " << endl;
624 }
956019b6 625
626 // It would be worth trying to calculate the inverse of the covariance matrix
627 // only once per fit, since it cannot change much in principle,
628 // and it would save a lot of computing time !!!!
04b5ea16 629
630 // Calculates Chi2
eb9c9dab 631 if (goodDeterminant) {
632 // with Multiple Scattering if inversion correct
d0bfce8d 633 // Inverse matrices without normalization
634 (*covBending) *= 1/normCovBending2;
635 (*covNonBending) *= 1/normCovNonBending2;
eb9c9dab 636 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
637 hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
638 hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
639 pbc1 = hit1->GetTrackParam()->GetBendingCoor();
640 hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
641 pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
642 for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
643 hit2 = (AliMUONTrackHit*)
644 (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
645 hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor();
646 pbc2 = hit2->GetTrackParam()->GetBendingCoor();
647 hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor();
648 pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor();
04b5ea16 649 Chi2 = Chi2 +
eb9c9dab 650 ((*covBending)(hitNumber2, hitNumber1) *
651 (hbc1 - pbc1) * (hbc2 - pbc2)) +
04b5ea16 652 ((*covNonBending)(hitNumber2, hitNumber1) *
eb9c9dab 653 (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
04b5ea16 654 }
655 }
eb9c9dab 656 } else {
657 // without Multiple Scattering if inversion impossible
658 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
659 hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
660 hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
661 pbc1 = hit1->GetTrackParam()->GetBendingCoor();
662 hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
663 pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
664 Chi2 = Chi2 +
665 ((hbc1 - pbc1) * (hbc1 - pbc1) /
666 hit1->GetHitForRecPtr()->GetBendingReso2()) +
667 ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) /
668 hit1->GetHitForRecPtr()->GetNonBendingReso2());
04b5ea16 669 }
670 }
671
672 delete covBending;
673 delete covNonBending;
eb9c9dab 674 delete [] msa2;
04b5ea16 675}
676
677Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
678{
679 // Returns square of multiple Coulomb scattering angle
680 // at TrackHit pointed to by "TrackHit"
681 Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
682 Double_t varMultipleScatteringAngle;
956019b6 683 AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
04b5ea16 684 AliMUONTrackParam *param = TrackHit->GetTrackParam();
956019b6 685 // Better implementation in AliMUONTrack class ????
04b5ea16 686 slopeBending = param->GetBendingSlope();
687 slopeNonBending = param->GetNonBendingSlope();
688 // thickness in radiation length for the current track,
689 // taking local angle into account
690 radiationLength =
691 trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
692 TMath::Sqrt(1.0 +
693 slopeBending * slopeBending + slopeNonBending * slopeNonBending);
694 inverseBendingMomentum2 =
695 param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
696 inverseTotalMomentum2 =
956019b6 697 inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
04b5ea16 698 (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
699 varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
956019b6 700 // The velocity is assumed to be 1 !!!!
04b5ea16 701 varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
702 varMultipleScatteringAngle * varMultipleScatteringAngle;
703 return varMultipleScatteringAngle;
704}