stdlib.h needed for non-Linux compilation
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.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.2  2000/06/15 07:58:48  morsch
19 Code from MUON-dev joined
20
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
49 #include <stdlib.h>
50
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     (&param1)->ExtrapToZ(zHit); // extrapolation
302     hit->SetTrackParam(&param1);
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() - (&param1)->GetNonBendingCoor();
309     Double_t dY =
310       hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
311     Chi2 =
312       Chi2 +
313       dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
314       dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
315   }
316 }