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