STEER/AliAODTrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 //     AOD track implementation of AliVParticle
  20 //     Author: Markus Oldenburg, CERN
  21 //     Markus.Oldenburg@cern.ch
  22 //-------------------------------------------------------------------------
  23
  24 #include "AliLog.h"
  25 #include "AliAODTrack.h"
  26
  27 ClassImp(AliAODTrack)
  28
  29 //______________________________________________________________________________
  30 AliAODTrack::AliAODTrack() :
  31   AliVParticle(),
  32   fChi2perNDF(-999.),
  33   fChi2MatchTrigger(0.),
  34   fFlags(0),
  35   fLabel(-999),
  36   fITSMuonClusterMap(0),
  37   fFilterMap(0),
  38   fID(-999),
  39   fCharge(-99),
  40   fType(kUndef),
  41   fCovMatrix(NULL),
  42   fDetPid(NULL),
  43   fProdVertex(NULL)
  44 {
  45   // default constructor
  46
  47   SetP();
  48   SetPosition((Float_t*)NULL);
  49   SetXYAtDCA(-999., -999.);
  50   SetPxPyPzAtDCA(-999., -999., -999.);
  51   SetPID((Float_t*)NULL);
  52 }
  53
  54 //______________________________________________________________________________
  55 AliAODTrack::AliAODTrack(Short_t id,
  56                          Int_t label,
  57                          Double_t p[3],
  58                          Bool_t cartesian,
  59                          Double_t x[3],
  60                          Bool_t isDCA,
  61                          Double_t covMatrix[21],
  62                          Short_t charge,
  63                          UChar_t itsClusMap,
  64                          Double_t pid[10],
  65                          AliAODVertex *prodVertex,
  66                          Bool_t usedForVtxFit,
  67                          Bool_t usedForPrimVtxFit,
  68                          AODTrk_t ttype,
  69                          UInt_t selectInfo,
  70                          Float_t chi2perNDF) :
  71   AliVParticle(),
  72   fChi2perNDF(chi2perNDF),
  73   fChi2MatchTrigger(0.),
  74   fFlags(0),
  75   fLabel(label),
  76   fITSMuonClusterMap(0),
  77   fFilterMap(selectInfo),
  78   fID(id),
  79   fCharge(charge),
  80   fType(ttype),
  81   fCovMatrix(NULL),
  82   fDetPid(NULL),
  83   fProdVertex(prodVertex)
  84 {
  85   // constructor
  86
  87   SetP(p, cartesian);
  88   SetPosition(x, isDCA);
  89   SetXYAtDCA(-999., -999.);
  90   SetPxPyPzAtDCA(-999., -999., -999.);
  91   SetUsedForVtxFit(usedForVtxFit);
  92   SetUsedForPrimVtxFit(usedForPrimVtxFit);
  93   if(covMatrix) SetCovMatrix(covMatrix);
  94   SetPID(pid);
  95   SetITSClusterMap(itsClusMap);
  96 }
  97
  98 //______________________________________________________________________________
  99 AliAODTrack::AliAODTrack(Short_t id,
 100                          Int_t label,
 101                          Float_t p[3],
 102                          Bool_t cartesian,
 103                          Float_t x[3],
 104                          Bool_t isDCA,
 105                          Float_t covMatrix[21],
 106                          Short_t charge,
 107                          UChar_t itsClusMap,
 108                          Float_t pid[10],
 109                          AliAODVertex *prodVertex,
 110                          Bool_t usedForVtxFit,
 111                          Bool_t usedForPrimVtxFit,
 112                          AODTrk_t ttype,
 113                          UInt_t selectInfo,
 114                          Float_t chi2perNDF) :
 115   AliVParticle(),
 116   fChi2perNDF(chi2perNDF),
 117   fChi2MatchTrigger(0.),
 118   fFlags(0),
 119   fLabel(label),
 120   fITSMuonClusterMap(0),
 121   fFilterMap(selectInfo),
 122   fID(id),
 123   fCharge(charge),
 124   fType(ttype),
 125   fCovMatrix(NULL),
 126   fDetPid(NULL),
 127   fProdVertex(prodVertex)
 128 {
 129   // constructor
 130
 131   SetP(p, cartesian);
 132   SetPosition(x, isDCA);
 133   SetXYAtDCA(-999., -999.);
 134   SetPxPyPzAtDCA(-999., -999., -999.);
 135   SetUsedForVtxFit(usedForVtxFit);
 136   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 137   if(covMatrix) SetCovMatrix(covMatrix);
 138   SetPID(pid);
 139   SetITSClusterMap(itsClusMap);
 140 }
 141
 142 //______________________________________________________________________________
 143 AliAODTrack::~AliAODTrack()
 144 {
 145   // destructor
 146   delete fCovMatrix;
 147 }
 148
 149
 150 //______________________________________________________________________________
 151 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
 152   AliVParticle(trk),
 153   fChi2perNDF(trk.fChi2perNDF),
 154   fChi2MatchTrigger(trk.fChi2MatchTrigger),
 155   fFlags(trk.fFlags),
 156   fLabel(trk.fLabel),
 157   fITSMuonClusterMap(trk.fITSMuonClusterMap),
 158   fFilterMap(trk.fFilterMap),
 159   fID(trk.fID),
 160   fCharge(trk.fCharge),
 161   fType(trk.fType),
 162   fCovMatrix(NULL),
 163   fDetPid(NULL),
 164   fProdVertex(trk.fProdVertex)
 165 {
 166   // Copy constructor
 167
 168   trk.GetP(fMomentum);
 169   trk.GetPosition(fPosition);
 170   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 171   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 172   SetUsedForVtxFit(trk.GetUsedForVtxFit());
 173   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 174   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 175   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 176   SetPID(trk.fPID);
 177 }
 178
 179 //______________________________________________________________________________
 180 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
 181 {
 182   // Assignment operator
 183   if(this!=&trk) {
 184
 185     AliVParticle::operator=(trk);
 186
 187     trk.GetP(fMomentum);
 188     trk.GetPosition(fPosition);
 189     trk.GetPID(fPID);
 190
 191     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 192     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 193
 194     fChi2perNDF = trk.fChi2perNDF;
 195     fChi2MatchTrigger = trk.fChi2MatchTrigger;
 196
 197     fFlags = trk.fFlags;
 198     fLabel = trk.fLabel;
 199
 200     fITSMuonClusterMap = trk.fITSMuonClusterMap;
 201     fFilterMap = trk.fFilterMap;
 202
 203     fID = trk.fID;
 204
 205     fCharge = trk.fCharge;
 206     fType = trk.fType;
 207
 208     delete fCovMatrix;
 209     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 210     else fCovMatrix=NULL;
 211     fProdVertex = trk.fProdVertex;
 212
 213     SetUsedForVtxFit(trk.GetUsedForVtxFit());
 214     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 215
 216     delete fDetPid;
 217     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 218     else fDetPid=NULL;
 219   }
 220
 221   return *this;
 222 }
 223
 224 //______________________________________________________________________________
 225 Double_t AliAODTrack::M(AODTrkPID_t pid) const
 226 {
 227   // Returns the mass.
 228   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
 229
 230   switch (pid) {
 231
 232   case kElectron :
 233     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
 234     break;
 235
 236   case kMuon :
 237     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
 238     break;
 239
 240   case kPion :
 241     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
 242     break;
 243
 244   case kKaon :
 245     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
 246     break;
 247
 248   case kProton :
 249     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
 250     break;
 251
 252   case kDeuteron :
 253     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
 254     break;
 255
 256   case kTriton :
 257     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
 258     break;
 259
 260   case kHelium3 :
 261     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
 262     break;
 263
 264   case kAlpha :
 265     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
 266     break;
 267
 268   case kUnknown :
 269     return -999.;
 270     break;
 271
 272   default :
 273     return -999.;
 274   }
 275 }
 276
 277 //______________________________________________________________________________
 278 Double_t AliAODTrack::E(AODTrkPID_t pid) const
 279 {
 280   // Returns the energy of the particle of a given pid.
 281
 282   if (pid != kUnknown) { // particle was identified
 283     Double_t m = M(pid);
 284     return TMath::Sqrt(P()*P() + m*m);
 285   } else { // pid unknown
 286     return -999.;
 287   }
 288 }
 289
 290 //______________________________________________________________________________
 291 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
 292 {
 293   // Returns the rapidity of a particle of a given pid.
 294
 295   if (pid != kUnknown) { // particle was identified
 296     Double_t e = E(pid);
 297     Double_t pz = Pz();
 298     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 299       return 0.5*TMath::Log((e+pz)/(e-pz));
 300     } else { // energy not known or equal to pz
 301       return -999.;
 302     }
 303   } else { // pid unknown
 304     return -999.;
 305   }
 306 }
 307
 308 //______________________________________________________________________________
 309 Double_t AliAODTrack::Y(Double_t m) const
 310 {
 311   // Returns the rapidity of a particle of a given mass.
 312
 313   if (m >= 0.) { // mass makes sense
 314     Double_t e = E(m);
 315     Double_t pz = Pz();
 316     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 317       return 0.5*TMath::Log((e+pz)/(e-pz));
 318     } else { // energy not known or equal to pz
 319       return -999.;
 320     }
 321   } else { // pid unknown
 322     return -999.;
 323   }
 324 }
 325
 326 //______________________________________________________________________________
 327 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const
 328 {
 329   // Returns the most probable PID array element.
 330
 331   Int_t nPID = 10;
 332   if (fPID) {
 333     AODTrkPID_t loc = kUnknown;
 334     Double_t max = 0.;
 335     Bool_t allTheSame = kTRUE;
 336
 337     for (Int_t iPID = 0; iPID < nPID; iPID++) {
 338       if (fPID[iPID] >= max) {
 339         if (fPID[iPID] > max) {
 340           allTheSame = kFALSE;
 341           max = fPID[iPID];
 342           loc = (AODTrkPID_t)iPID;
 343         } else {
 344           allTheSame = kTRUE;
 345         }
 346       }
 347     }
 348
 349     return allTheSame ? kUnknown : loc;
 350   } else {
 351     return kUnknown;
 352   }
 353 }
 354
 355 //______________________________________________________________________________
 356 void AliAODTrack::ConvertAliPIDtoAODPID()
 357 {
 358   // Converts AliPID array.
 359   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
 360   // Everything else has to be set to zero.
 361
 362   fPID[kDeuteron] = 0.;
 363   fPID[kTriton]   = 0.;
 364   fPID[kHelium3]  = 0.;
 365   fPID[kAlpha]    = 0.;
 366   fPID[kUnknown]  = 0.;
 367
 368   return;
 369 }
 370
 371
 372 //______________________________________________________________________________
 373 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian)
 374 {
 375   // Set the momentum
 376
 377   if (p) {
 378     if (cartesian) {
 379       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
 380       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
 381
 382       fMomentum[0] = TMath::Sqrt(pt2); // pt
 383       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
 384       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
 385     } else {
 386       fMomentum[0] = p[0];  // pt
 387       fMomentum[1] = p[1];  // phi
 388       fMomentum[2] = p[2];  // theta
 389     }
 390   } else {
 391     fMomentum[0] = -999.;
 392     fMomentum[1] = -999.;
 393     fMomentum[2] = -999.;
 394   }
 395 }
 396
 397 //______________________________________________________________________________
 398 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca)
 399 {
 400   // set the position
 401
 402   if (x) {
 403     if (!dca) {
 404       ResetBit(kIsDCA);
 405
 406       fPosition[0] = x[0];
 407       fPosition[1] = x[1];
 408       fPosition[2] = x[2];
 409     } else {
 410       SetBit(kIsDCA);
 411       // don't know any better yet
 412       fPosition[0] = -999.;
 413       fPosition[1] = -999.;
 414       fPosition[2] = -999.;
 415     }
 416   } else {
 417     ResetBit(kIsDCA);
 418
 419     fPosition[0] = -999.;
 420     fPosition[1] = -999.;
 421     fPosition[2] = -999.;
 422   }
 423 }
 424
 425 //______________________________________________________________________________
 426 void AliAODTrack::SetDCA(Double_t d, Double_t z)
 427 {
 428   // set the dca
 429   fPosition[0] = d;
 430   fPosition[1] = z;
 431   fPosition[2] = 0.;
 432   SetBit(kIsDCA);
 433 }
 434
 435 //______________________________________________________________________________
 436 void AliAODTrack::Print(Option_t* /* option */) const
 437 {
 438   // prints information about AliAODTrack
 439
 440   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle());
 441   printf("        px = %f\n", Px());
 442   printf("        py = %f\n", Py());
 443   printf("        pz = %f\n", Pz());
 444   printf("        pt = %f\n", Pt());
 445   printf("      1/pt = %f\n", OneOverPt());
 446   printf("     theta = %f\n", Theta());
 447   printf("       phi = %f\n", Phi());
 448   printf("  chi2/NDF = %f\n", Chi2perNDF());
 449   printf("    charge = %d\n", Charge());
 450 }
 451
 452 void AliAODTrack::SetMatchTrigger(Int_t matchTrig){
 453 //
 454 // Set the MUON trigger information
 455   switch(matchTrig){
 456     case 0: // 0 track does not match trigger
 457       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 458       break;
 459     case 1: // 1 track match but does not pass pt cut
 460       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
 461       break;
 462     case 2: // 2 track match Low pt cut
 463       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
 464       break;
 465     case 3: // 3 track match High pt cut
 466       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
 467       break;
 468     default:
 469       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 470       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
 471   }
 472 }
 473
 474 Int_t AliAODTrack::HitsMT(Int_t istation, Int_t iplane, Char_t *cathode){
 475 //
 476 // Retrieve hit information for MUON identified by  (station, plane, cathode)
 477   if(cathode){
 478     if(cathode[0]=='x'||cathode[0]=='X'){
 479       if(istation==1){
 480         if(iplane==1)
 481           return (fITSMuonClusterMap&0x8000)?1:0;
 482         else if(iplane==2)
 483           return (fITSMuonClusterMap&0x4000)?1:0;
 484         else
 485           return 0;
 486       }else if(istation==2){
 487         if(iplane==1)
 488           return (fITSMuonClusterMap&0x2000)?1:0;
 489         else if(iplane==2)
 490           return (fITSMuonClusterMap&0x1000)?1:0;
 491         else
 492           return 0;
 493       }else{
 494         return 0;
 495       }
 496     }else if(cathode[0]=='y'||cathode[0]=='Y'){
 497       if(istation==1){
 498         if(iplane==1)
 499           return (fITSMuonClusterMap&0x0800)?1:0;
 500         else if(iplane==2)
 501           return (fITSMuonClusterMap&0x0400)?1:0;
 502         else
 503           return 0;
 504       }else if(istation==2){
 505         if(iplane==1)
 506           return (fITSMuonClusterMap&0x0200)?1:0;
 507         else if(iplane==2)
 508           return (fITSMuonClusterMap&0x0100)?1:0;
 509         else
 510           return 0;
 511       }else{
 512         return 0;
 513       }
 514     }else{
 515       return 0;
 516     }
 517   }else{
 518     if(istation==1){
 519       if(iplane==1)
 520         return (HitsMT(1,1,"X")||HitsMT(1,1,"Y"))?1:0;
 521       else if(iplane==2)
 522         return (HitsMT(1,2,"X")||HitsMT(1,2,"Y"))?1:0;
 523       else
 524         return 0;
 525     }else if(istation==2){
 526       if(iplane==1)
 527         return (HitsMT(2,1,"X")||HitsMT(2,1,"Y"))?1:0;
 528       else if(iplane==2)
 529         return (HitsMT(2,2,"X")||HitsMT(2,2,"Y"))?1:0;
 530       else
 531         return 0;
 532     }else{
 533       return 0;
 534     }
 535   }
 536 }
 537
 538 Int_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber){
 539   //
 540   // Retrieve hit information for MUON Tracker/Trigger Chamber
 541   // WARNING: chamber number start from 1 instead of 0
 542
 543   if (MuonChamber > 0 && MuonChamber < 11) {
 544     return ((GetMUONClusterMap() & BIT(MuonChamber-1)) != 0) ? 1 : 0;
 545   } else {
 546     switch(MuonChamber){
 547       case 11:
 548         return HitsMT(1,1);
 549       case 12:
 550         return HitsMT(1,2);
 551       case 13:
 552         return HitsMT(2,1);
 553       case 14:
 554         return HitsMT(2,2);
 555       default:
 556         printf("Unknown MUON chamber: %d\n",MuonChamber);
 557         return 0;
 558     }
 559   }
 560 }
 561
 562 Bool_t AliAODTrack::PropagateTo(Double_t xk, Double_t b) {
 563   //----------------------------------------------------------------
 564   // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
 565   // This is in local coordinates!!!
 566   //----------------------------------------------------------------
 567
 568   Double_t alpha = 0.;
 569   Double_t localP[3] = {Px(), Py(), Pz()}; // set global (sic!) p
 570   Global2LocalMomentum(localP, Charge(), alpha); // convert global to local momentum
 571
 572   AliAODVertex *origin = (AliAODVertex*)fProdVertex.GetObject();
 573   Double_t localX[3] = {origin->GetX(), origin->GetY(), origin->GetZ()}; // set global (sic!) location of first track point
 574   Global2LocalPosition(localX, alpha); // convert global to local position
 575
 576   Double_t &fX = localX[0];
 577
 578   Double_t dx=xk-fX;
 579   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
 580
 581   Double_t crv=localP[0]*b*kB2C;
 582   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
 583
 584   Double_t f1=localP[1], f2=f1 + crv*dx;
 585   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
 586   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
 587
 588   Double_t &fP0=localX[1], &fP1=localX[2], &fP2=localP[0], &fP3=localP[1], &fP4=localP[2];
 589   /* covariance matrix to be fixed!
 590   Double_t
 591   &fC00=fC[0],
 592   &fC10=fC[1],   &fC11=fC[2],
 593   &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
 594   &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],
 595   &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];
 596   */
 597   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
 598
 599   fX=xk;
 600   fP0 += dx*(f1+f2)/(r1+r2);
 601   fP1 += dx*(r2 + f2*(f1+f2)/(r1+r2))*fP3;
 602   fP2 += dx*crv;
 603
 604   //f = F - 1
 605
 606   //Double_t f02=    dx/(r1*r1*r1);
 607   Double_t cc=crv/fP4;
 608   Double_t f04=0.5*dx*dx/(r1*r1*r1);         f04*=cc;
 609   //Double_t f12=    dx*fP3*f1/(r1*r1*r1);
 610   Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);  f14*=cc;
 611   //Double_t f13=    dx/r1;
 612   Double_t f24=    dx;                       f24*=cc;
 613
 614   /* covariance matrix to be fixed!
 615   //b = C*ft
 616   Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
 617   Double_t b02=f24*fC40;
 618   Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
 619   Double_t b12=f24*fC41;
 620   Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
 621   Double_t b22=f24*fC42;
 622   Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
 623   Double_t b42=f24*fC44;
 624   Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
 625   Double_t b32=f24*fC43;
 626
 627   //a = f*b = f*C*ft
 628   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
 629   Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
 630   Double_t a22=f24*b42;
 631
 632   //F*C*Ft = C + (b + bt + a)
 633   fC00 += b00 + b00 + a00;
 634   fC10 += b10 + b01 + a01;
 635   fC20 += b20 + b02 + a02;
 636   fC30 += b30;
 637   fC40 += b40;
 638   fC11 += b11 + b11 + a11;
 639   fC21 += b21 + b12 + a12;
 640   fC31 += b31;
 641   fC41 += b41;
 642   fC22 += b22 + b22 + a22;
 643   fC32 += b32;
 644   fC42 += b42;
 645   */
 646
 647   Local2GlobalMomentum(localP, alpha); // convert local to global momentum
 648   SetP(localP);
 649
 650   return kTRUE;
 651 }
 652