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 }