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