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