1 #include "AliAODPair.h"
2 //_________________________________________________________________________
3 ///////////////////////////////////////////////////////////////////////////
7 // class implements pair of particles and taking care of caluclation (almost)
8 // all of pair properties (Qinv, InvMass,...)
10 // more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
12 ////////////////////////////////////////////////////////////////////////////
14 #include "AliVAODParticle.h"
15 #include "AliTrackPoints.h"
18 /************************************************************************/
19 AliAODPair::AliAODPair(Bool_t rev):
24 fQSideLCMSNotCalc(kTRUE),
26 fQOutLCMSNotCalc(kTRUE),
28 fQLongLCMSNotCalc(kTRUE),
30 fQtLCMSNotCalc(kTRUE),
36 fInvMassNotCalc(kTRUE),
44 fKStarCompNotCalc(kTRUE),
52 fMassSqrNotCalc(kTRUE),
55 fAvarageDistance(0.0),
56 fAvarageDistanceNotCalc(kTRUE),
68 fGammaLCMSNotCalc(kTRUE),
71 //value of rev defines if it is Swapped
72 //if you pass kTRUE swpaped pair will NOT be created
73 //though you wont be able to get the swaped pair from this pair
75 if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
78 /************************************************************************/
80 AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
85 fQSideLCMSNotCalc(kTRUE),
87 fQOutLCMSNotCalc(kTRUE),
89 fQLongLCMSNotCalc(kTRUE),
91 fQtLCMSNotCalc(kTRUE),
97 fInvMassNotCalc(kTRUE),
101 fKStarNotCalc(kTRUE),
105 fKStarCompNotCalc(kTRUE),
113 fMassSqrNotCalc(kTRUE),
115 fQInvLNotCalc(kTRUE),
116 fAvarageDistance(0.0),
117 fAvarageDistanceNotCalc(kTRUE),
127 fDiffsNotCalc(kTRUE),
129 fGammaLCMSNotCalc(kTRUE),
132 //value of rev defines if it is Swapped
133 //if you pass kTRUE swpaped pair will NOT be created
134 //though you wont be able to get the swaped pair from this pair
136 if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
139 /************************************************************************/
140 AliAODPair::AliAODPair(const AliAODPair& in):
146 fQSideLCMSNotCalc(kTRUE),
148 fQOutLCMSNotCalc(kTRUE),
150 fQLongLCMSNotCalc(kTRUE),
152 fQtLCMSNotCalc(kTRUE),
158 fInvMassNotCalc(kTRUE),
162 fKStarNotCalc(kTRUE),
166 fKStarCompNotCalc(kTRUE),
174 fMassSqrNotCalc(kTRUE),
176 fQInvLNotCalc(kTRUE),
177 fAvarageDistance(0.0),
178 fAvarageDistanceNotCalc(kTRUE),
188 fDiffsNotCalc(kTRUE),
190 fGammaLCMSNotCalc(kTRUE),
196 /************************************************************************/
198 AliAODPair& AliAODPair::operator=(const AliAODPair& in)
204 /************************************************************************/
206 Double_t AliAODPair::GetInvMass()
208 //Returns qinv value for a pair
211 CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
213 if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
214 else fInvMass = TMath::Sqrt(fInvMassSqr);
216 fInvMassNotCalc = kFALSE;
220 /************************************************************************/
222 Double_t AliAODPair::GetQSideLCMS()
224 //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
226 if (fQSideLCMSNotCalc)
228 fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
229 fQSideLCMSNotCalc = kFALSE;
233 /************************************************************************/
235 Double_t AliAODPair::GetQOutLCMS()
237 //caculates Qout in Center Of Mass Longitudionally Co-Moving
243 if (fPart1->Mass() != fPart2->Mass())
247 Double_t beta = fPzSum/fESum;
248 Double_t gamma = GetGammaToLCMS();
249 Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
250 Double_t x = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
251 beta = 2.0*GetKt()/GetMt();
252 gamma = GetMt()/GetQInv();
253 fQOutLCMS = gamma * (x - beta * el);
256 //beta=fPzSum/fESum; // Longit. V == beta
257 Double_t beta=fPzSum/fESum;
258 Double_t gamma = GetGammaToLCMS();
260 Double_t cosphi=fPxSum/(2.0*GetKt()); // cos(phi)
261 Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
263 // ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
264 // ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
266 tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
267 Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
268 Double_t part1Px = tmp;
270 tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
271 Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
272 Double_t part2Px = tmp;
275 // LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
276 // LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
277 Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
278 Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
280 Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
281 Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
282 Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
283 Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
284 Double_t sumE=part1E+part2E;
285 Double_t sumPx=part1Px+part2Px;
286 Double_t sumPy=part1Py+part2Py;
287 Double_t sumPZ=part1Pz+part2Pz;
288 Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
290 Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
291 Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
292 fQOutLCMS=(part1Px-part2Px);//== id
293 fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
297 Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
298 fQOutLCMS = 0.5*k2/GetKt();
299 // if (non-id) fQOutLCMS=fQOutLCMS - sumPx*HF;
303 fQOutLCMSNotCalc = kFALSE;
307 /************************************************************************/
309 Double_t AliAODPair::GetQLongLCMS()
311 //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
312 if (fQLongLCMSNotCalc)
316 Double_t beta = fPzSum/fESum;
317 fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
318 fQLongLCMSNotCalc = kFALSE;
322 /************************************************************************/
324 Double_t AliAODPair::GetQtLCMS()
326 //returns Q transverse CMS longitudionally co-moving
329 fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
330 fQtLCMSNotCalc = kFALSE;
334 /************************************************************************/
336 Double_t AliAODPair::GetQt()
338 //returns Q transverse CMS longitudionally co-moving
344 Double_t dotprod = fPxSum*fPxDiff + fPySum*fPyDiff + fPzSum*fPzDiff;
345 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
346 klen = TMath::Sqrt(klen);
347 Double_t qlen = fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff;
348 qlen = TMath::Sqrt(qlen);
349 Double_t tmp = klen*qlen;
356 Double_t cosopenangle = dotprod/tmp;
357 Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
359 fQt = sinopenangle*qlen;
364 /************************************************************************/
366 Double_t AliAODPair::GetKt()
368 //calculates the evarage momentum of the pair
372 fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
377 /************************************************************************/
379 Double_t AliAODPair::GetKStar()
381 //calculates invariant velocity difference
386 Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
387 Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
390 Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
391 Error("GetKStar","Particle1:");
393 Error("GetKStar","Particle2:");
395 Error("GetKStar","");
398 fKStarNotCalc = kFALSE;
401 Double_t pinv = TMath::Sqrt(mtrans - ptrans);
403 Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
410 Info("GetKStar","Sqrt of negative number q = %f",q);
411 Error("GetKStar","Particle1:");
413 Error("GetKStar","Particle2:");
416 fKStarNotCalc = kFALSE;
422 fKStarNotCalc = kFALSE;
426 /************************************************************************/
427 Double_t AliAODPair::GetKStarOut()
429 CalculateKStarComp();
432 /************************************************************************/
433 Double_t AliAODPair::GetKStarSide()
435 CalculateKStarComp();
438 /************************************************************************/
439 Double_t AliAODPair::GetKStarLong()
441 CalculateKStarComp();
444 /************************************************************************/
446 Double_t AliAODPair::GetQInv()
449 //warning for non-id particles you want to use 2*KStar
453 fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
454 fQInvNotCalc = kFALSE;
458 /************************************************************************/
460 Double_t AliAODPair::GetGammaToLCMS()
462 //calculates gamma factor of the boost to LCMS
463 if(fGammaLCMSNotCalc)
466 Double_t beta = fPzSum/fESum;
467 fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
468 fGammaLCMSNotCalc = kFALSE;
472 /************************************************************************/
474 Double_t AliAODPair::GetGammaToTransverse()
476 //calculates gamma factor of the boost to LCMS
477 Double_t beta = 2.0*GetKt() / GetMt();
478 Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
482 /************************************************************************/
484 Double_t AliAODPair::GetMt()
486 //Calculates transverse mass of the pair
490 fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
495 /************************************************************************/
497 Double_t AliAODPair::GetAvarageDistance()
499 //returns and buffers avarage distance between two tracks calculated
500 // out of track points (see AliAODTrackPoints class)
502 if (fAvarageDistanceNotCalc)
504 fAvarageDistance = AvDistance();
505 fAvarageDistanceNotCalc = kFALSE;
507 return fAvarageDistance;
509 /************************************************************************/
511 Double_t AliAODPair::AvDistance()
513 //returns avarage distance between two tracks in range
514 //as defined in Track-Points of AliVAODParticle
515 //returns negative value if error uccured f.g. tracks do not have track-points
516 AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
518 {//it could be simulated pair
519 // Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
523 AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
526 // Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
530 return tpts1->AvarageDistance(*tpts2);
532 /************************************************************************/
534 Double_t AliAODPair::GetR()
536 //Returns distance between particles vertexes in thir CMS
540 Double_t vxDiff = fPart1->Vx() - fPart2->Vx();
541 Double_t vyDiff = fPart1->Vy() - fPart2->Vy();
542 Double_t vzDiff = fPart1->Vz() - fPart2->Vz();
544 Double_t r = TMath::Sqrt( vxDiff*vxDiff + vyDiff*vyDiff + vzDiff*vzDiff );
548 /************************************************************************/
550 Double_t AliAODPair::GetRStar()
552 //Returns distance between particles vertexes in thir CMS
557 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
558 klen = TMath::Sqrt(klen);
560 Double_t aBeta = klen/fESum;
561 Double_t aGamma = 1.0/TMath::Sqrt(1.0 - aBeta*aBeta);
564 Double_t alpha = -TMath::ATan2(fPySum,fPzSum);
565 Double_t beta = TMath::ATan2(fPxSum,TMath::Hypot(fPySum,fPzSum));
567 Double_t sinalpha = TMath::Sin(alpha);
568 Double_t cosalpha = TMath::Cos(alpha);
569 Double_t sinbeta = TMath::Sin(beta);
570 Double_t cosbeta = TMath::Cos(beta);
572 Double_t v1xP = fPart1->Vx();
573 Double_t v2xP = fPart2->Vx();
574 Double_t v1yP = fPart1->Vy()*cosalpha + fPart1->Vz()*sinalpha;
575 Double_t v2yP = fPart2->Vy()*cosalpha + fPart2->Vz()*sinalpha;
576 Double_t v1zP =-fPart1->Vy()*sinalpha + fPart1->Vz()*cosalpha;
577 Double_t v2zP =-fPart2->Vy()*sinalpha + fPart2->Vz()*cosalpha;
580 ///////////////////////////////////////////////////
582 // Double_t p1yP = fPart1->Py()*cosalpha + fPart1->Pz()*sinalpha;
583 // Double_t p2yP = fPart2->Py()*cosalpha + fPart2->Pz()*sinalpha;
585 // Double_t p1zP =-fPart1->Py()*sinalpha + fPart1->Pz()*cosalpha;
586 // Double_t p2zP =-fPart2->Py()*sinalpha + fPart2->Pz()*cosalpha;
589 // Double_t p1x = fPart1->Px()*cosbeta - p1zP*sinbeta;
590 // Double_t p2x = fPart2->Px()*cosbeta - p2zP*sinbeta;
591 // Double_t p1z = fPart1->Px()*sinbeta + p1zP*cosbeta;
592 // Double_t p2z = fPart2->Px()*sinbeta + p2zP*cosbeta;
594 // Info("","%f %f %f",p1yP,p2yP,p1yP+p2yP);
595 // Info("","%f %f %f",p1x,p2x,p1x+p2x);
597 // Info("","%f %f ",p1x+p2x,p1yP+p2yP);
599 ///////////////////////////////////////////////////
602 Double_t v1x = v1xP*cosbeta - v1zP*sinbeta;
603 Double_t v2x = v2xP*cosbeta - v2zP*sinbeta;
606 Double_t v1z = v1xP*sinbeta + v1zP*cosbeta;
607 Double_t v2z = v2xP*sinbeta + v2zP*cosbeta;
610 Double_t v1zB=aGamma*(v1z-aBeta*fPart1->T());
611 Double_t v2zB=aGamma*(v2z-aBeta*fPart2->T());
615 Double_t dx = v1x - v2x;
616 Double_t dy = v1y - v2y;
617 Double_t dz = v1zB - v2zB;
619 Double_t rstar = TMath::Sqrt( dx*dx + dy*dy + dz*dz);
623 /************************************************************************/
625 void AliAODPair::MirrorSecond()
627 //makes local copy of the second particle and mirrors their momenta
628 //for its deletion is responsible who calls this method
629 fPart2 = (AliVAODParticle*)fPart2->Clone();
630 fPart2->SetMomentum(-fPart2->Px(),-fPart2->Py(),-fPart2->Pz(),fPart2->E());
633 /************************************************************************/
635 void AliAODPair::DeleteSecond()
637 //Deletes second particle
642 void AliAODPair::Print()
644 if (fPart1) fPart1->Print();
645 if (fPart2) fPart2->Print();
647 Info("Print","GetKStar() %f",GetKStar());
648 Info("Print","GetKt() %f",GetKt() );
649 Info("Print","QInv %f", GetQInv() );
650 Info("Print","GetQOutLCMS() %f",GetQOutLCMS() );
651 Info("Print","GetQSideLCMS %f",GetQSideLCMS() );
652 Info("Print","GetQLongLCMS() %f", GetQLongLCMS());
653 Info("Print","GetDeltaTheta() %f", GetDeltaTheta());
654 Info("Print","GetDeltaPhi() %f", GetDeltaPhi());