| | 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 | /* $Id$ */ |
| | 17 | |
| | 18 | /////////////////////////////////////////////////// |
| | 19 | // |
| | 20 | // Reconstructed track |
| | 21 | // in |
| | 22 | // ALICE |
| | 23 | // dimuon |
| | 24 | // spectrometer |
| | 25 | // |
| | 26 | /////////////////////////////////////////////////// |
| | 27 | |
| | 28 | #include "AliMUONTrack.h" |
| | 29 | |
| | 30 | #include "AliMUONTrackReconstructor.h" |
| | 31 | #include "AliMUONHitForRec.h" |
| | 32 | #include "AliMUONSegment.h" |
| | 33 | #include "AliMUONTrackHit.h" |
| | 34 | #include "AliMUONTriggerTrack.h" |
| | 35 | #include "AliMUONConstants.h" |
| | 36 | |
| | 37 | #include "AliLog.h" |
| | 38 | |
| | 39 | #include <Riostream.h> // for cout |
| | 40 | #include <TMath.h> |
| | 41 | #include <TMatrixD.h> |
| | 42 | #include <TObjArray.h> |
| | 43 | #include <TVirtualFitter.h> |
| | 44 | |
| | 45 | #include <stdlib.h> // for exit() |
| | 46 | |
| | 47 | // Functions to be minimized with Minuit |
| | 48 | void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); |
| | 49 | void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); |
| | 50 | |
| | 51 | void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); |
| | 52 | |
| | 53 | Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit); |
| | 54 | |
| | 55 | TVirtualFitter* AliMUONTrack::fgFitter = NULL; |
| | 56 | |
| | 57 | ClassImp(AliMUONTrack) // Class implementation in ROOT context |
| | 58 | |
| | 59 | //__________________________________________________________________________ |
| | 60 | AliMUONTrack::AliMUONTrack() |
| | 61 | : TObject() |
| | 62 | { |
| | 63 | // Default constructor |
| | 64 | fgFitter = 0; |
| | 65 | fTrackReconstructor = 0; |
| | 66 | fTrackHitsPtr = NULL; |
| | 67 | fTrackParamAtHit = NULL; |
| | 68 | fHitForRecAtHit = NULL; |
| | 69 | fTrackID = 0; |
| | 70 | } |
| | 71 | |
| | 72 | //__________________________________________________________________________ |
| | 73 | AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONTrackReconstructor* TrackReconstructor) |
| | 74 | : TObject() |
| | 75 | { |
| | 76 | // Constructor from two Segment's |
| | 77 | fTrackReconstructor = TrackReconstructor; // link back to TrackReconstructor |
| | 78 | // memory allocation for the TObjArray of pointers to reconstructed TrackHit's |
| | 79 | fTrackHitsPtr = new TObjArray(10); |
| | 80 | fNTrackHits = 0; |
| | 81 | if (BegSegment) { //AZ |
| | 82 | AddSegment(BegSegment); // add hits from BegSegment |
| | 83 | AddSegment(EndSegment); // add hits from EndSegment |
| | 84 | fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z |
| | 85 | SetTrackParamAtVertex(); // set track parameters at vertex |
| | 86 | } |
| | 87 | fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); |
| | 88 | fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); |
| | 89 | // set fit conditions... |
| | 90 | fFitMCS = 0; |
| | 91 | fFitNParam = 3; |
| | 92 | fFitStart = 1; |
| | 93 | fFitFMin = -1.0; |
| | 94 | fMatchTrigger = kFALSE; |
| | 95 | fChi2MatchTrigger = 0; |
| | 96 | fTrackID = 0; |
| | 97 | return; |
| | 98 | } |
| | 99 | |
| | 100 | //__________________________________________________________________________ |
| | 101 | AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONTrackReconstructor* TrackReconstructor) |
| | 102 | : TObject() |
| | 103 | { |
| | 104 | // Constructor from one Segment and one HitForRec |
| | 105 | fTrackReconstructor = TrackReconstructor; // link back to TrackReconstructor |
| | 106 | // memory allocation for the TObjArray of pointers to reconstructed TrackHit's |
| | 107 | fTrackHitsPtr = new TObjArray(10); |
| | 108 | fNTrackHits = 0; |
| | 109 | AddSegment(Segment); // add hits from Segment |
| | 110 | AddHitForRec(HitForRec); // add HitForRec |
| | 111 | fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z |
| | 112 | SetTrackParamAtVertex(); // set track parameters at vertex |
| | 113 | fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); |
| | 114 | fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); |
| | 115 | // set fit conditions... |
| | 116 | fFitMCS = 0; |
| | 117 | fFitNParam = 3; |
| | 118 | fFitStart = 1; |
| | 119 | fFitFMin = -1.0; |
| | 120 | fMatchTrigger = kFALSE; |
| | 121 | fChi2MatchTrigger = 0; |
| | 122 | fTrackID = 0; |
| | 123 | return; |
| | 124 | } |
| | 125 | |
| | 126 | //__________________________________________________________________________ |
| | 127 | AliMUONTrack::~AliMUONTrack() |
| | 128 | { |
| | 129 | // Destructor |
| | 130 | if (fTrackHitsPtr) { |
| | 131 | // delete the TObjArray of pointers to TrackHit's |
| | 132 | delete fTrackHitsPtr; |
| | 133 | fTrackHitsPtr = NULL; |
| | 134 | } |
| | 135 | |
| | 136 | if (fTrackParamAtHit) { |
| | 137 | // delete the TClonesArray of pointers to TrackParam |
| | 138 | delete fTrackParamAtHit; |
| | 139 | fTrackParamAtHit = NULL; |
| | 140 | } |
| | 141 | |
| | 142 | if (fHitForRecAtHit) { |
| | 143 | // delete the TClonesArray of pointers to HitForRec |
| | 144 | delete fHitForRecAtHit; |
| | 145 | fHitForRecAtHit = NULL; |
| | 146 | } |
| | 147 | } |
| | 148 | |
| | 149 | //__________________________________________________________________________ |
| | 150 | AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack) |
| | 151 | : TObject(theMUONTrack) |
| | 152 | { |
| | 153 | //fTrackReconstructor = new AliMUONTrackReconstructor(*MUONTrack.fTrackReconstructor); |
| | 154 | // is it right ? |
| | 155 | // NO, because it would use dummy copy constructor |
| | 156 | // and AliMUONTrack is not the owner of its TrackReconstructor |
| | 157 | fTrackReconstructor = theMUONTrack.fTrackReconstructor; |
| | 158 | fTrackParamAtVertex = theMUONTrack.fTrackParamAtVertex; |
| | 159 | |
| | 160 | // necessary to make a copy of the objects and not only the pointers in TObjArray. |
| | 161 | fTrackHitsPtr = new TObjArray(10); |
| | 162 | for (Int_t index = 0; index < (theMUONTrack.fTrackHitsPtr)->GetEntriesFast(); index++) { |
| | 163 | AliMUONTrackHit *trackHit = new AliMUONTrackHit(*(AliMUONTrackHit*)(theMUONTrack.fTrackHitsPtr)->At(index)); |
| | 164 | fTrackHitsPtr->Add(trackHit); |
| | 165 | } |
| | 166 | fTrackHitsPtr->SetOwner(); // nedeed for deleting TClonesArray |
| | 167 | |
| | 168 | // necessary to make a copy of the objects and not only the pointers in TClonesArray. |
| | 169 | fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); |
| | 170 | for (Int_t index = 0; index < (theMUONTrack.fTrackParamAtHit)->GetEntriesFast(); index++) { |
| | 171 | {new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) |
| | 172 | AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));} |
| | 173 | } |
| | 174 | |
| | 175 | // necessary to make a copy of the objects and not only the pointers in TClonesArray. |
| | 176 | fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); |
| | 177 | for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) { |
| | 178 | {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) |
| | 179 | AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));} |
| | 180 | } |
| | 181 | |
| | 182 | fNTrackHits = theMUONTrack.fNTrackHits; |
| | 183 | fFitMCS = theMUONTrack.fFitMCS; |
| | 184 | fFitNParam = theMUONTrack.fFitNParam; |
| | 185 | fFitFMin = theMUONTrack.fFitFMin; |
| | 186 | fFitStart = theMUONTrack.fFitStart; |
| | 187 | fMatchTrigger = theMUONTrack.fMatchTrigger; |
| | 188 | fChi2MatchTrigger = theMUONTrack.fChi2MatchTrigger; |
| | 189 | fTrackID = theMUONTrack.fTrackID; |
| | 190 | } |
| | 191 | |
| | 192 | //__________________________________________________________________________ |
| | 193 | AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& theMUONTrack) |
| | 194 | { |
| | 195 | |
| | 196 | // check assignement to self |
| | 197 | if (this == &theMUONTrack) |
| | 198 | return *this; |
| | 199 | |
| | 200 | // base class assignement |
| | 201 | TObject::operator=(theMUONTrack); |
| | 202 | |
| | 203 | // fTrackReconstructor = new AliMUONTrackReconstructor(*MUONTrack.fTrackReconstructor); |
| | 204 | // is it right ? |
| | 205 | // is it right ? NO because it would use dummy copy constructor |
| | 206 | fTrackReconstructor = theMUONTrack.fTrackReconstructor; |
| | 207 | fTrackParamAtVertex = theMUONTrack.fTrackParamAtVertex; |
| | 208 | |
| | 209 | // necessary to make a copy of the objects and not only the pointers in TObjArray. |
| | 210 | fTrackHitsPtr = new TObjArray(10); |
| | 211 | for (Int_t index = 0; index < (theMUONTrack.fTrackHitsPtr)->GetEntriesFast(); index++) { |
| | 212 | AliMUONTrackHit *trackHit = new AliMUONTrackHit(*(AliMUONTrackHit*)(theMUONTrack.fTrackHitsPtr)->At(index)); |
| | 213 | fTrackHitsPtr->Add(trackHit); |
| | 214 | } |
| | 215 | fTrackHitsPtr->SetOwner(); // nedeed for deleting TClonesArray |
| | 216 | |
| | 217 | // necessary to make a copy of the objects and not only the pointers in TClonesArray. |
| | 218 | fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); |
| | 219 | for (Int_t index = 0; index < (theMUONTrack.fTrackParamAtHit)->GetEntriesFast(); index++) { |
| | 220 | {new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) |
| | 221 | AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));} |
| | 222 | } |
| | 223 | |
| | 224 | // necessary to make a copy of the objects and not only the pointers in TClonesArray. |
| | 225 | fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); |
| | 226 | for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) { |
| | 227 | {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) |
| | 228 | AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));} |
| | 229 | } |
| | 230 | |
| | 231 | fNTrackHits = theMUONTrack.fNTrackHits; |
| | 232 | fFitMCS = theMUONTrack.fFitMCS; |
| | 233 | fFitNParam = theMUONTrack.fFitNParam; |
| | 234 | fFitFMin = theMUONTrack.fFitFMin; |
| | 235 | fFitStart = theMUONTrack.fFitStart; |
| | 236 | fMatchTrigger = theMUONTrack.fMatchTrigger; |
| | 237 | fChi2MatchTrigger = theMUONTrack.fChi2MatchTrigger; |
| | 238 | fTrackID = theMUONTrack.fTrackID; |
| | 239 | |
| | 240 | return *this; |
| | 241 | } |
| | 242 | |
| | 243 | //__________________________________________________________________________ |
| | 244 | void AliMUONTrack::Remove() |
| | 245 | { |
| | 246 | // Remove current track from array of tracks, |
| | 247 | // and corresponding track hits from array of track hits. |
| | 248 | // Compress the TClonesArray it belongs to. |
| | 249 | AliMUONTrackHit *nextTrackHit; |
| | 250 | AliMUONTrackReconstructor *eventRec = this->fTrackReconstructor; |
| | 251 | TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr(); |
| | 252 | // Loop over all track hits of track |
| | 253 | AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First(); |
| | 254 | while (trackHit) { |
| | 255 | nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit); |
| | 256 | // Remove TrackHit from event TClonesArray. |
| | 257 | // Destructor is called, |
| | 258 | // hence links between HitForRec's and TrackHit's are updated |
| | 259 | trackHitsPtr->Remove(trackHit); |
| | 260 | trackHit = nextTrackHit; |
| | 261 | } |
| | 262 | // Remove the track from event TClonesArray |
| | 263 | // Destructor is called, |
| | 264 | // hence space for TObjArray of pointers to TrackHit's is freed |
| | 265 | eventRec->GetRecTracksPtr()->Remove(this); |
| | 266 | // Number of tracks decreased by 1 |
| | 267 | eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1); |
| | 268 | // Compress event TClonesArray of Track's: |
| | 269 | // this is essential to retrieve the TClonesArray afterwards |
| | 270 | eventRec->GetRecTracksPtr()->Compress(); |
| | 271 | // Compress event TClonesArray of TrackHit's: |
| | 272 | // this is probably also essential to retrieve the TClonesArray afterwards |
| | 273 | trackHitsPtr->Compress(); |
| | 274 | } |
| | 275 | |
| | 276 | //__________________________________________________________________________ |
| | 277 | void AliMUONTrack::SetFitMCS(Int_t FitMCS) |
| | 278 | { |
| | 279 | // Set multiple Coulomb scattering option for track fit "fFitMCS" |
| | 280 | // from "FitMCS" argument: 0 without, 1 with |
| | 281 | if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS; |
| | 282 | // better implementation with enum(with, without) ???? |
| | 283 | else { |
| | 284 | cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl; |
| | 285 | cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl; |
| | 286 | exit(0); |
| | 287 | } |
| | 288 | return; |
| | 289 | } |
| | 290 | |
| | 291 | //__________________________________________________________________________ |
| | 292 | void AliMUONTrack::SetFitNParam(Int_t FitNParam) |
| | 293 | { |
| | 294 | // Set number of parameters for track fit "fFitNParam" from "FitNParam": |
| | 295 | // 3 for momentum, 5 for momentum and position |
| | 296 | if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam; |
| | 297 | else { |
| | 298 | cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl; |
| | 299 | cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl; |
| | 300 | exit(0); |
| | 301 | } |
| | 302 | return; |
| | 303 | } |
| | 304 | |
| | 305 | //__________________________________________________________________________ |
| | 306 | void AliMUONTrack::SetFitStart(Int_t FitStart) |
| | 307 | { |
| | 308 | // Set multiple Coulomb scattering option for track fit "fFitStart" |
| | 309 | // from "FitStart" argument: 0 without, 1 with |
| | 310 | if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart; |
| | 311 | // better implementation with enum(vertex, firstHit) ???? |
| | 312 | else { |
| | 313 | cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl; |
| | 314 | cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl; |
| | 315 | exit(0); |
| | 316 | } |
| | 317 | return; |
| | 318 | } |
| | 319 | |
| | 320 | //__________________________________________________________________________ |
| | 321 | AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const { |
| | 322 | // Get pointer to TrackParamAtFirstHit |
| | 323 | return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} |
| | 324 | |
| | 325 | //__________________________________________________________________________ |
| | 326 | void AliMUONTrack::RecursiveDump(void) const |
| | 327 | { |
| | 328 | // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's |
| | 329 | AliMUONTrackHit *trackHit; |
| | 330 | AliMUONHitForRec *hitForRec; |
| | 331 | cout << "Recursive dump of Track: " << this << endl; |
| | 332 | // Track |
| | 333 | this->Dump(); |
| | 334 | for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) { |
| | 335 | trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]); |
| | 336 | // TrackHit |
| | 337 | cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl; |
| | 338 | trackHit->Dump(); |
| | 339 | hitForRec = trackHit->GetHitForRecPtr(); |
| | 340 | // HitForRec |
| | 341 | cout << "HitForRec: " << hitForRec << endl; |
| | 342 | hitForRec->Dump(); |
| | 343 | } |
| | 344 | return; |
| | 345 | } |
| | 346 | |
| | 347 | //__________________________________________________________________________ |
| | 348 | Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut) const |
| | 349 | { |
| | 350 | // Return kTRUE/kFALSE for each chamber if hit is compatible or not |
| | 351 | TClonesArray *hitArray, *thisHitArray; |
| | 352 | AliMUONHitForRec *hit, *thisHit; |
| | 353 | Int_t chamberNumber; |
| | 354 | Float_t deltaZ; |
| | 355 | Float_t deltaZMax = 1.; // 1 cm |
| | 356 | Float_t chi2 = 0; |
| | 357 | Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()]; |
| | 358 | |
| | 359 | for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { |
| | 360 | nCompHit[ch] = kFALSE; |
| | 361 | } |
| | 362 | |
| | 363 | thisHitArray = this->GetHitForRecAtHit(); |
| | 364 | |
| | 365 | hitArray = Track->GetHitForRecAtHit(); |
| | 366 | |
| | 367 | for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) { |
| | 368 | thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis); |
| | 369 | chamberNumber = thisHit->GetChamberNumber(); |
| | 370 | if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue; |
| | 371 | nCompHit[chamberNumber] = kFALSE; |
| | 372 | for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) { |
| | 373 | hit = (AliMUONHitForRec*) hitArray->At(iH); |
| | 374 | deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ()); |
| | 375 | chi2 = thisHit->NormalizedChi2WithHitForRec(hit,Sigma2Cut); // set cut to 4 sigmas |
| | 376 | if (chi2 < 3. && deltaZ < deltaZMax) { |
| | 377 | nCompHit[chamberNumber] = kTRUE; |
| | 378 | break; |
| | 379 | } |
| | 380 | } |
| | 381 | } |
| | 382 | |
| | 383 | return nCompHit; |
| | 384 | } |
| | 385 | |
| | 386 | //__________________________________________________________________________ |
| | 387 | Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const |
| | 388 | { |
| | 389 | // Returns the number of hits in common |
| | 390 | // between the current track ("this") |
| | 391 | // and the track pointed to by "Track". |
| | 392 | Int_t hitsInCommon = 0; |
| | 393 | AliMUONTrackHit *trackHit1, *trackHit2; |
| | 394 | // Loop over hits of first track |
| | 395 | trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First(); |
| | 396 | while (trackHit1) { |
| | 397 | // Loop over hits of second track |
| | 398 | trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First(); |
| | 399 | while (trackHit2) { |
| | 400 | // Increment "hitsInCommon" if both TrackHits point to the same HitForRec |
| | 401 | if ( (trackHit1->GetHitForRecPtr()) == |
| | 402 | (trackHit2->GetHitForRecPtr()) ) hitsInCommon++; |
| | 403 | trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2); |
| | 404 | } // trackHit2 |
| | 405 | trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1); |
| | 406 | } // trackHit1 |
| | 407 | return hitsInCommon; |
| | 408 | } |
| | 409 | |
| | 410 | //__________________________________________________________________________ |
| | 411 | void AliMUONTrack::MatchTriggerTrack(TClonesArray *triggerTrackArray) |
| | 412 | { |
| | 413 | // Match this track with one trigger track if possible |
| | 414 | AliMUONTrackParam trackParam; |
| | 415 | AliMUONTriggerTrack *triggerTrack; |
| | 416 | Double_t xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2; |
| | 417 | Double_t nSigmaCut2; |
| | 418 | |
| | 419 | Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY |
| | 420 | Double_t distTriggerTrack[3] = {0,0,0}; |
| | 421 | |
| | 422 | fMatchTrigger = kFALSE; |
| | 423 | fChi2MatchTrigger = 0; |
| | 424 | |
| | 425 | trackParam = *((AliMUONTrackParam*) fTrackParamAtHit->Last()); |
| | 426 | trackParam.ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber |
| | 427 | |
| | 428 | nSigmaCut2 = fTrackReconstructor->GetMaxSigma2Distance(); // nb of sigma**2 for cut |
| | 429 | xTrack = trackParam.GetNonBendingCoor(); |
| | 430 | yTrack = trackParam.GetBendingCoor(); |
| | 431 | ySlopeTrack = trackParam.GetBendingSlope(); |
| | 432 | dTrigTrackMin2 = 999; |
| | 433 | |
| | 434 | triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->First(); |
| | 435 | while(triggerTrack){ |
| | 436 | distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0]; |
| | 437 | distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1]; |
| | 438 | distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2]; |
| | 439 | dTrigTrack2 = 0; |
| | 440 | for (Int_t iVar = 0; iVar < 3; iVar++) |
| | 441 | dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar]; |
| | 442 | if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < nSigmaCut2) { |
| | 443 | dTrigTrackMin2 = dTrigTrack2; |
| | 444 | fMatchTrigger = kTRUE; |
| | 445 | fChi2MatchTrigger = dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY) |
| | 446 | } |
| | 447 | triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->After(triggerTrack); |
| | 448 | } |
| | 449 | |
| | 450 | } |
| | 451 | //__________________________________________________________________________ |
| | 452 | void AliMUONTrack::Fit() |
| | 453 | { |
| | 454 | // Fit the current track ("this"), |
| | 455 | // with or without multiple Coulomb scattering according to "fFitMCS", |
| | 456 | // with the number of parameters given by "fFitNParam" |
| | 457 | // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary), |
| | 458 | // starting, according to "fFitStart", |
| | 459 | // with track parameters at vertex or at the first TrackHit. |
| | 460 | // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before |
| | 461 | // by calling the corresponding Set methods. |
| | 462 | Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y; |
| | 463 | char parName[50]; |
| | 464 | AliMUONTrackParam *trackParam; |
| | 465 | // Check if Minuit is initialized... |
| | 466 | fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ??? |
| | 467 | fgFitter->Clear(); // necessary ???? probably yes |
| | 468 | // how to reset the printout number at every fit ???? |
| | 469 | // is there any risk to leave it like that ???? |
| | 470 | // how to go faster ???? choice of Minuit parameters like EDM ???? |
| | 471 | // choice of function to be minimized according to fFitMCS |
| | 472 | if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2); |
| | 473 | else fgFitter->SetFCN(TrackChi2MCS); |
| | 474 | // Switch off printout |
| | 475 | arg[0] = -1; |
| | 476 | fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!! |
| | 477 | // No warnings |
| | 478 | fgFitter->ExecuteCommand("SET NOW", arg, 0); |
| | 479 | // Parameters according to "fFitStart" |
| | 480 | // (should be a function to be used at every place where needed ????) |
| | 481 | if (fFitStart == 0) trackParam = &fTrackParamAtVertex; |
| | 482 | else trackParam = this->GetTrackParamAtFirstHit(); |
| | 483 | // set first 3 Minuit parameters |
| | 484 | // could be tried with no limits for the search (min=max=0) ???? |
| | 485 | fgFitter->SetParameter(0, "InvBenP", |
| | 486 | trackParam->GetInverseBendingMomentum(), |
| | 487 | 0.003, -0.4, 0.4); |
| | 488 | fgFitter->SetParameter(1, "BenS", |
| | 489 | trackParam->GetBendingSlope(), |
| | 490 | 0.001, -0.5, 0.5); |
| | 491 | fgFitter->SetParameter(2, "NonBenS", |
| | 492 | trackParam->GetNonBendingSlope(), |
| | 493 | 0.001, -0.5, 0.5); |
| | 494 | if (fFitNParam == 5) { |
| | 495 | // set last 2 Minuit parameters |
| | 496 | // mandatory limits in Bending to avoid NaN values of parameters |
| | 497 | fgFitter->SetParameter(3, "X", |
| | 498 | trackParam->GetNonBendingCoor(), |
| | 499 | 0.03, -500.0, 500.0); |
| | 500 | // mandatory limits in non Bending to avoid NaN values of parameters |
| | 501 | fgFitter->SetParameter(4, "Y", |
| | 502 | trackParam->GetBendingCoor(), |
| | 503 | 0.10, -500.0, 500.0); |
| | 504 | } |
| | 505 | // search without gradient calculation in the function |
| | 506 | fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0); |
| | 507 | // minimization |
| | 508 | fgFitter->ExecuteCommand("MINIMIZE", arg, 0); |
| | 509 | // exit from Minuit |
| | 510 | // fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ???? |
| | 511 | // get results into "invBenP", "benC", "nonBenC" ("x", "y") |
| | 512 | fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper); |
| | 513 | fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper); |
| | 514 | fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper); |
| | 515 | if (fFitNParam == 5) { |
| | 516 | fgFitter->GetParameter(3, parName, x, errorParam, lower, upper); |
| | 517 | fgFitter->GetParameter(4, parName, y, errorParam, lower, upper); |
| | 518 | } |
| | 519 | // result of the fit into track parameters |
| | 520 | trackParam->SetInverseBendingMomentum(invBenP); |
| | 521 | trackParam->SetBendingSlope(benC); |
| | 522 | trackParam->SetNonBendingSlope(nonBenC); |
| | 523 | if (fFitNParam == 5) { |
| | 524 | trackParam->SetNonBendingCoor(x); |
| | 525 | trackParam->SetBendingCoor(y); |
| | 526 | } |
| | 527 | // global result of the fit |
| | 528 | Double_t fedm, errdef; |
| | 529 | Int_t npari, nparx; |
| | 530 | fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx); |
| | 531 | } |
| | 532 | |
| | 533 | //__________________________________________________________________________ |
| | 534 | void AliMUONTrack::AddSegment(AliMUONSegment* Segment) |
| | 535 | { |
| | 536 | // Add Segment to the track |
| | 537 | AddHitForRec(Segment->GetHitForRec1()); // 1st hit |
| | 538 | AddHitForRec(Segment->GetHitForRec2()); // 2nd hit |
| | 539 | } |
| | 540 | |
| | 541 | //__________________________________________________________________________ |
| | 542 | void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) |
| | 543 | { |
| | 544 | // Add HitForRec to the track: |
| | 545 | // actual TrackHit into TClonesArray of TrackHit's for the event; |
| | 546 | // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track |
| | 547 | TClonesArray *recTrackHitsPtr = this->fTrackReconstructor->GetRecTrackHitsPtr(); |
| | 548 | Int_t eventTrackHits = this->fTrackReconstructor->GetNRecTrackHits(); |
| | 549 | // event |
| | 550 | AliMUONTrackHit* trackHit = |
| | 551 | new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec); |
| | 552 | this->fTrackReconstructor->SetNRecTrackHits(eventTrackHits + 1); |
| | 553 | // track |
| | 554 | if (fTrackHitsPtr->IsOwner()) AliFatal("fTrackHitsPtr is owner"); |
| | 555 | fTrackHitsPtr->Add(trackHit); |
| | 556 | fNTrackHits++; |
| | 557 | } |
| | 558 | |
| | 559 | //__________________________________________________________________________ |
| | 560 | void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const |
| | 561 | { |
| | 562 | // Set track parameters at TrackHit with index "indexHit" |
| | 563 | // from the track parameters pointed to by "TrackParam". |
| | 564 | //PH AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); |
| | 565 | AliMUONTrackHit* trackHit = (AliMUONTrackHit*) (fTrackHitsPtr->At(indexHit)); |
| | 566 | trackHit->SetTrackParam(TrackParam); |
| | 567 | } |
| | 568 | |
| | 569 | //__________________________________________________________________________ |
| | 570 | void AliMUONTrack::SetTrackParamAtVertex() |
| | 571 | { |
| | 572 | // Set track parameters at vertex. |
| | 573 | // TrackHit's are assumed to be only in stations(1..) 4 and 5, |
| | 574 | // and sorted according to increasing Z.. |
| | 575 | // Parameters are calculated from information in HitForRec's |
| | 576 | // of first and last TrackHit's. |
| | 577 | AliMUONTrackParam *trackParam = |
| | 578 | &fTrackParamAtVertex; // pointer to track parameters |
| | 579 | // Pointer to HitForRec of first TrackHit |
| | 580 | AliMUONHitForRec *firstHit = |
| | 581 | ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr(); |
| | 582 | // Pointer to HitForRec of last TrackHit |
| | 583 | AliMUONHitForRec *lastHit = |
| | 584 | ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr(); |
| | 585 | // Z difference between first and last hits |
| | 586 | Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ(); |
| | 587 | // bending slope in stations(1..) 4 and 5 |
| | 588 | Double_t bendingSlope = |
| | 589 | (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ; |
| | 590 | trackParam->SetBendingSlope(bendingSlope); |
| | 591 | // impact parameter |
| | 592 | Double_t impactParam = |
| | 593 | firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and lastHit ???? |
| | 594 | // signed bending momentum |
| | 595 | Double_t signedBendingMomentum = |
| | 596 | fTrackReconstructor->GetBendingMomentumFromImpactParam(impactParam); |
| | 597 | trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum); |
| | 598 | // bending slope at vertex |
| | 599 | trackParam-> |
| | 600 | SetBendingSlope(bendingSlope + |
| | 601 | impactParam / fTrackReconstructor->GetSimpleBPosition()); |
| | 602 | // non bending slope |
| | 603 | Double_t nonBendingSlope = |
| | 604 | (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ; |
| | 605 | trackParam->SetNonBendingSlope(nonBendingSlope); |
| | 606 | // vertex coordinates at (0,0,0) |
| | 607 | trackParam->SetZ(0.0); |
| | 608 | trackParam->SetBendingCoor(0.0); |
| | 609 | trackParam->SetNonBendingCoor(0.0); |
| | 610 | } |
| | 611 | |
| | 612 | //__________________________________________________________________________ |
| | 613 | void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) |
| | 614 | { |
| | 615 | // Return the "Chi2" to be minimized with Minuit for track fitting, |
| | 616 | // with "NParam" parameters |
| | 617 | // and their current values in array pointed to by "Param". |
| | 618 | // Assumes that the track hits are sorted according to increasing Z. |
| | 619 | // Track parameters at each TrackHit are updated accordingly. |
| | 620 | // Multiple Coulomb scattering is not taken into account |
| | 621 | AliMUONTrack *trackBeingFitted; |
| | 622 | AliMUONTrackHit* hit; |
| | 623 | AliMUONTrackParam param1; |
| | 624 | Int_t hitNumber; |
| | 625 | Double_t zHit; |
| | 626 | Chi2 = 0.0; // initialize Chi2 |
| | 627 | // copy of track parameters to be fitted |
| | 628 | trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); |
| | 629 | if (trackBeingFitted->GetFitStart() == 0) |
| | 630 | param1 = *(trackBeingFitted->GetTrackParamAtVertex()); |
| | 631 | else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); |
| | 632 | // Minuit parameters to be fitted into this copy |
| | 633 | param1.SetInverseBendingMomentum(Param[0]); |
| | 634 | param1.SetBendingSlope(Param[1]); |
| | 635 | param1.SetNonBendingSlope(Param[2]); |
| | 636 | if (NParam == 5) { |
| | 637 | param1.SetNonBendingCoor(Param[3]); |
| | 638 | param1.SetBendingCoor(Param[4]); |
| | 639 | } |
| | 640 | // Follow track through all planes of track hits |
| | 641 | for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) { |
| | 642 | hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; |
| | 643 | zHit = hit->GetHitForRecPtr()->GetZ(); |
| | 644 | // do something special if 2 hits with same Z ???? |
| | 645 | // security against infinite loop ???? |
| | 646 | (¶m1)->ExtrapToZ(zHit); // extrapolation |
| | 647 | hit->SetTrackParam(¶m1); |
| | 648 | // Increment Chi2 |
| | 649 | // done hit per hit, with hit resolution, |
| | 650 | // and not with point and angle like in "reco_muon.F" !!!! |
| | 651 | // Needs to add multiple scattering contribution ???? |
| | 652 | Double_t dX = |
| | 653 | hit->GetHitForRecPtr()->GetNonBendingCoor() - (¶m1)->GetNonBendingCoor(); |
| | 654 | Double_t dY = |
| | 655 | hit->GetHitForRecPtr()->GetBendingCoor() - (¶m1)->GetBendingCoor(); |
| | 656 | Chi2 = |
| | 657 | Chi2 + |
| | 658 | dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() + |
| | 659 | dY * dY / hit->GetHitForRecPtr()->GetBendingReso2(); |
| | 660 | } |
| | 661 | } |
| | 662 | |
| | 663 | //__________________________________________________________________________ |
| | 664 | void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) |
| | 665 | { |
| | 666 | // Return the "Chi2" to be minimized with Minuit for track fitting, |
| | 667 | // with "NParam" parameters |
| | 668 | // and their current values in array pointed to by "Param". |
| | 669 | // Assumes that the track hits are sorted according to increasing Z. |
| | 670 | // Track parameters at each TrackHit are updated accordingly. |
| | 671 | // Multiple Coulomb scattering is taken into account with covariance matrix. |
| | 672 | AliMUONTrack *trackBeingFitted; |
| | 673 | AliMUONTrackParam param1; |
| | 674 | Chi2 = 0.0; // initialize Chi2 |
| | 675 | // copy of track parameters to be fitted |
| | 676 | trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); |
| | 677 | if (trackBeingFitted->GetFitStart() == 0) |
| | 678 | param1 = *(trackBeingFitted->GetTrackParamAtVertex()); |
| | 679 | else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); |
| | 680 | // Minuit parameters to be fitted into this copy |
| | 681 | param1.SetInverseBendingMomentum(Param[0]); |
| | 682 | param1.SetBendingSlope(Param[1]); |
| | 683 | param1.SetNonBendingSlope(Param[2]); |
| | 684 | if (NParam == 5) { |
| | 685 | param1.SetNonBendingCoor(Param[3]); |
| | 686 | param1.SetBendingCoor(Param[4]); |
| | 687 | } |
| | 688 | |
| | 689 | AliMUONTrackHit *hit; |
| | 690 | Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; |
| | 691 | Double_t z, z1, z2, z3; |
| | 692 | AliMUONTrackHit *hit1, *hit2, *hit3; |
| | 693 | Double_t hbc1, hbc2, pbc1, pbc2; |
| | 694 | Double_t hnbc1, hnbc2, pnbc1, pnbc2; |
| | 695 | Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); |
| | 696 | TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit); |
| | 697 | TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit); |
| | 698 | Double_t *msa2 = new Double_t[numberOfHit]; |
| | 699 | |
| | 700 | // Predicted coordinates and multiple scattering angles are first calculated |
| | 701 | for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { |
| | 702 | hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; |
| | 703 | z = hit->GetHitForRecPtr()->GetZ(); |
| | 704 | // do something special if 2 hits with same Z ???? |
| | 705 | // security against infinite loop ???? |
| | 706 | (¶m1)->ExtrapToZ(z); // extrapolation |
| | 707 | hit->SetTrackParam(¶m1); |
| | 708 | // square of multiple scattering angle at current hit, with one chamber |
| | 709 | msa2[hitNumber] = MultipleScatteringAngle2(hit); |
| | 710 | // correction for eventual missing hits or multiple hits in a chamber, |
| | 711 | // according to the number of chambers |
| | 712 | // between the current hit and the previous one |
| | 713 | chCurrent = hit->GetHitForRecPtr()->GetChamberNumber(); |
| | 714 | if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev); |
| | 715 | chPrev = chCurrent; |
| | 716 | } |
| | 717 | |
| | 718 | // Calculates the covariance matrix |
| | 719 | for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { |
| | 720 | hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; |
| | 721 | z1 = hit1->GetHitForRecPtr()->GetZ(); |
| | 722 | for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { |
| | 723 | hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; |
| | 724 | z2 = hit2->GetHitForRecPtr()->GetZ(); |
| | 725 | // initialization to 0 (diagonal plus upper triangular part) |
| | 726 | (*covBending)(hitNumber2, hitNumber1) = 0.0; |
| | 727 | // contribution from multiple scattering in bending plane: |
| | 728 | // loop over upstream hits |
| | 729 | for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) { |
| | 730 | hit3 = (AliMUONTrackHit*) |
| | 731 | (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3]; |
| | 732 | z3 = hit3->GetHitForRecPtr()->GetZ(); |
| | 733 | (*covBending)(hitNumber2, hitNumber1) = |
| | 734 | (*covBending)(hitNumber2, hitNumber1) + |
| | 735 | ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); |
| | 736 | } |
| | 737 | // equal contribution from multiple scattering in non bending plane |
| | 738 | (*covNonBending)(hitNumber2, hitNumber1) = |
| | 739 | (*covBending)(hitNumber2, hitNumber1); |
| | 740 | if (hitNumber1 == hitNumber2) { |
| | 741 | // Diagonal elements: add contribution from position measurements |
| | 742 | // in bending plane |
| | 743 | (*covBending)(hitNumber2, hitNumber1) = |
| | 744 | (*covBending)(hitNumber2, hitNumber1) + |
| | 745 | hit1->GetHitForRecPtr()->GetBendingReso2(); |
| | 746 | // and in non bending plane |
| | 747 | (*covNonBending)(hitNumber2, hitNumber1) = |
| | 748 | (*covNonBending)(hitNumber2, hitNumber1) + |
| | 749 | hit1->GetHitForRecPtr()->GetNonBendingReso2(); |
| | 750 | } |
| | 751 | else { |
| | 752 | // Non diagonal elements: symmetrization |
| | 753 | // for bending plane |
| | 754 | (*covBending)(hitNumber1, hitNumber2) = |
| | 755 | (*covBending)(hitNumber2, hitNumber1); |
| | 756 | // and non bending plane |
| | 757 | (*covNonBending)(hitNumber1, hitNumber2) = |
| | 758 | (*covNonBending)(hitNumber2, hitNumber1); |
| | 759 | } |
| | 760 | } // for (hitNumber2 = hitNumber1;... |
| | 761 | } // for (hitNumber1 = 0;... |
| | 762 | |
| | 763 | // Inversion of covariance matrices |
| | 764 | // with "mnvertLocal", local "mnvert" function of Minuit. |
| | 765 | // One cannot use directly "mnvert" since "TVirtualFitter" does not know it. |
| | 766 | // One will have to replace this local function by the right inversion function |
| | 767 | // from a specialized Root package for symmetric positive definite matrices, |
| | 768 | // when available!!!! |
| | 769 | Int_t ifailBending; |
| | 770 | mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, |
| | 771 | ifailBending); |
| | 772 | Int_t ifailNonBending; |
| | 773 | mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, |
| | 774 | ifailNonBending); |
| | 775 | |
| | 776 | // It would be worth trying to calculate the inverse of the covariance matrix |
| | 777 | // only once per fit, since it cannot change much in principle, |
| | 778 | // and it would save a lot of computing time !!!! |
| | 779 | |
| | 780 | // Calculates Chi2 |
| | 781 | if ((ifailBending == 0) && (ifailNonBending == 0)) { |
| | 782 | // with Multiple Scattering if inversion correct |
| | 783 | for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { |
| | 784 | hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; |
| | 785 | hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); |
| | 786 | pbc1 = hit1->GetTrackParam()->GetBendingCoor(); |
| | 787 | hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); |
| | 788 | pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); |
| | 789 | for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) { |
| | 790 | hit2 = (AliMUONTrackHit*) |
| | 791 | (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; |
| | 792 | hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor(); |
| | 793 | pbc2 = hit2->GetTrackParam()->GetBendingCoor(); |
| | 794 | hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor(); |
| | 795 | pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor(); |
| | 796 | Chi2 = Chi2 + |
| | 797 | ((*covBending)(hitNumber2, hitNumber1) * |
| | 798 | (hbc1 - pbc1) * (hbc2 - pbc2)) + |
| | 799 | ((*covNonBending)(hitNumber2, hitNumber1) * |
| | 800 | (hnbc1 - pnbc1) * (hnbc2 - pnbc2)); |
| | 801 | } |
| | 802 | } |
| | 803 | } else { |
| | 804 | // without Multiple Scattering if inversion impossible |
| | 805 | for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { |
| | 806 | hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; |
| | 807 | hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); |
| | 808 | pbc1 = hit1->GetTrackParam()->GetBendingCoor(); |
| | 809 | hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); |
| | 810 | pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); |
| | 811 | Chi2 = Chi2 + |
| | 812 | ((hbc1 - pbc1) * (hbc1 - pbc1) / |
| | 813 | hit1->GetHitForRecPtr()->GetBendingReso2()) + |
| | 814 | ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / |
| | 815 | hit1->GetHitForRecPtr()->GetNonBendingReso2()); |
| | 816 | } |
| | 817 | } |
| | 818 | |
| | 819 | delete covBending; |
| | 820 | delete covNonBending; |
| | 821 | delete [] msa2; |
| | 822 | } |
| | 823 | |
| | 824 | Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) |
| | 825 | { |
| | 826 | // Returns square of multiple Coulomb scattering angle |
| | 827 | // at TrackHit pointed to by "TrackHit" |
| | 828 | Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; |
| | 829 | Double_t varMultipleScatteringAngle; |
| | 830 | AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); |
| | 831 | AliMUONTrackParam *param = TrackHit->GetTrackParam(); |
| | 832 | // Better implementation in AliMUONTrack class ???? |
| | 833 | slopeBending = param->GetBendingSlope(); |
| | 834 | slopeNonBending = param->GetNonBendingSlope(); |
| | 835 | // thickness in radiation length for the current track, |
| | 836 | // taking local angle into account |
| | 837 | radiationLength = |
| | 838 | trackBeingFitted->GetTrackReconstructor()->GetChamberThicknessInX0() * |
| | 839 | TMath::Sqrt(1.0 + |
| | 840 | slopeBending * slopeBending + slopeNonBending * slopeNonBending); |
| | 841 | inverseBendingMomentum2 = |
| | 842 | param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum(); |
| | 843 | inverseTotalMomentum2 = |
| | 844 | inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) / |
| | 845 | (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); |
| | 846 | varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength)); |
| | 847 | // The velocity is assumed to be 1 !!!! |
| | 848 | varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * |
| | 849 | varMultipleScatteringAngle * varMultipleScatteringAngle; |
| | 850 | return varMultipleScatteringAngle; |
| | 851 | } |
| | 852 | |
| | 853 | //______________________________________________________________________________ |
| | 854 | void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail) |
| | 855 | { |
| | 856 | //*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-* |
| | 857 | //*-* ========================== |
| | 858 | //*-* inverts a symmetric matrix. matrix is first scaled to |
| | 859 | //*-* have all ones on the diagonal (equivalent to change of units) |
| | 860 | //*-* but no pivoting is done since matrix is positive-definite. |
| | 861 | //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* |
| | 862 | |
| | 863 | // taken from TMinuit package of Root (l>=n) |
| | 864 | // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp |
| | 865 | // Double_t localVERTs[n], localVERTq[n], localVERTpp[n]; |
| | 866 | Double_t * localVERTs = new Double_t[n]; |
| | 867 | Double_t * localVERTq = new Double_t[n]; |
| | 868 | Double_t * localVERTpp = new Double_t[n]; |
| | 869 | // fMaxint changed to localMaxint |
| | 870 | Int_t localMaxint = n; |
| | 871 | |
| | 872 | /* System generated locals */ |
| | 873 | Int_t aOffset; |
| | 874 | |
| | 875 | /* Local variables */ |
| | 876 | Double_t si; |
| | 877 | Int_t i, j, k, kp1, km1; |
| | 878 | |
| | 879 | /* Parameter adjustments */ |
| | 880 | aOffset = l + 1; |
| | 881 | a -= aOffset; |
| | 882 | |
| | 883 | /* Function Body */ |
| | 884 | ifail = 0; |
| | 885 | if (n < 1) goto L100; |
| | 886 | if (n > localMaxint) goto L100; |
| | 887 | //*-*- scale matrix by sqrt of diag elements |
| | 888 | for (i = 1; i <= n; ++i) { |
| | 889 | si = a[i + i*l]; |
| | 890 | if (si <= 0) goto L100; |
| | 891 | localVERTs[i-1] = 1 / TMath::Sqrt(si); |
| | 892 | } |
| | 893 | for (i = 1; i <= n; ++i) { |
| | 894 | for (j = 1; j <= n; ++j) { |
| | 895 | a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1]; |
| | 896 | } |
| | 897 | } |
| | 898 | //*-*- . . . start main loop . . . . |
| | 899 | for (i = 1; i <= n; ++i) { |
| | 900 | k = i; |
| | 901 | //*-*- preparation for elimination step1 |
| | 902 | if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l]; |
| | 903 | else goto L100; |
| | 904 | localVERTpp[k-1] = 1; |
| | 905 | a[k + k*l] = 0; |
| | 906 | kp1 = k + 1; |
| | 907 | km1 = k - 1; |
| | 908 | if (km1 < 0) goto L100; |
| | 909 | else if (km1 == 0) goto L50; |
| | 910 | else goto L40; |
| | 911 | L40: |
| | 912 | for (j = 1; j <= km1; ++j) { |
| | 913 | localVERTpp[j-1] = a[j + k*l]; |
| | 914 | localVERTq[j-1] = a[j + k*l]*localVERTq[k-1]; |
| | 915 | a[j + k*l] = 0; |
| | 916 | } |
| | 917 | L50: |
| | 918 | if (k - n < 0) goto L51; |
| | 919 | else if (k - n == 0) goto L60; |
| | 920 | else goto L100; |
| | 921 | L51: |
| | 922 | for (j = kp1; j <= n; ++j) { |
| | 923 | localVERTpp[j-1] = a[k + j*l]; |
| | 924 | localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1]; |
| | 925 | a[k + j*l] = 0; |
| | 926 | } |
| | 927 | //*-*- elimination proper |
| | 928 | L60: |
| | 929 | for (j = 1; j <= n; ++j) { |
| | 930 | for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; } |
| | 931 | } |
| | 932 | } |
| | 933 | //*-*- elements of left diagonal and unscaling |
| | 934 | for (j = 1; j <= n; ++j) { |
| | 935 | for (k = 1; k <= j; ++k) { |
| | 936 | a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1]; |
| | 937 | a[j + k*l] = a[k + j*l]; |
| | 938 | } |
| | 939 | } |
| | 940 | delete [] localVERTs; |
| | 941 | delete [] localVERTq; |
| | 942 | delete [] localVERTpp; |
| | 943 | return; |
| | 944 | //*-*- failure return |
| | 945 | L100: |
| | 946 | delete [] localVERTs; |
| | 947 | delete [] localVERTq; |
| | 948 | delete [] localVERTpp; |
| | 949 | ifail = 1; |
| | 950 | } /* mnvertLocal */ |
| | 951 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental