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),
95 fInvMassNotCalc(kTRUE),
103 fKStarCompNotCalc(kTRUE),
111 fMassSqrNotCalc(kTRUE),
113 fQInvLNotCalc(kTRUE),
114 fAvarageDistance(0.0),
115 fAvarageDistanceNotCalc(kTRUE),
125 fDiffsNotCalc(kTRUE),
127 fGammaLCMSNotCalc(kTRUE),
130 //value of rev defines if it is Swapped
131 //if you pass kTRUE swpaped pair will NOT be created
132 //though you wont be able to get the swaped pair from this pair
134 if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
137 /************************************************************************/
138 AliAODPair::AliAODPair(const AliAODPair& in):
144 fQSideLCMSNotCalc(kTRUE),
146 fQOutLCMSNotCalc(kTRUE),
148 fQLongLCMSNotCalc(kTRUE),
150 fQtLCMSNotCalc(kTRUE),
154 fInvMassNotCalc(kTRUE),
158 fKStarNotCalc(kTRUE),
162 fKStarCompNotCalc(kTRUE),
170 fMassSqrNotCalc(kTRUE),
172 fQInvLNotCalc(kTRUE),
173 fAvarageDistance(0.0),
174 fAvarageDistanceNotCalc(kTRUE),
184 fDiffsNotCalc(kTRUE),
186 fGammaLCMSNotCalc(kTRUE),
192 /************************************************************************/
194 AliAODPair& AliAODPair::operator=(const AliAODPair& in)
200 /************************************************************************/
202 Double_t AliAODPair::GetInvMass()
204 //Returns qinv value for a pair
207 CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
209 if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
210 else fInvMass = TMath::Sqrt(fInvMassSqr);
212 fInvMassNotCalc = kFALSE;
216 /************************************************************************/
218 Double_t AliAODPair::GetQSideLCMS()
220 //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
222 if (fQSideLCMSNotCalc)
224 fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
225 fQSideLCMSNotCalc = kFALSE;
229 /************************************************************************/
231 Double_t AliAODPair::GetQOutLCMS()
233 //caculates Qout in Center Of Mass Longitudionally Co-Moving
239 if (fPart1->Mass() != fPart2->Mass())
243 Double_t beta = fPzSum/fESum;
244 Double_t gamma = GetGammaToLCMS();
245 Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
246 Double_t x = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
247 beta = 2.0*GetKt()/GetMt();
248 gamma = GetMt()/GetQInv();
249 fQOutLCMS = gamma * (x - beta * el);
252 //beta=fPzSum/fESum; // Longit. V == beta
253 Double_t beta=fPzSum/fESum;
254 Double_t gamma = GetGammaToLCMS();
256 Double_t cosphi=fPxSum/(2.0*GetKt()); // cos(phi)
257 Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
259 // ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
260 // ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
262 tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
263 Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
264 Double_t part1Px = tmp;
266 tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
267 Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
268 Double_t part2Px = tmp;
271 // LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
272 // LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
273 Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
274 Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
276 Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
277 Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
278 Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
279 Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
280 Double_t sumE=part1E+part2E;
281 Double_t sumPx=part1Px+part2Px;
282 Double_t sumPy=part1Py+part2Py;
283 Double_t sumPZ=part1Pz+part2Pz;
284 Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
286 Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
287 Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
288 fQOutLCMS=(part1Px-part2Px);//== id
289 fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
293 Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
294 fQOutLCMS = 0.5*k2/GetKt();
295 // if (non-id) fQOutLCMS=fQOutLCMS - sumPx*HF;
299 fQOutLCMSNotCalc = kFALSE;
303 /************************************************************************/
305 Double_t AliAODPair::GetQLongLCMS()
307 //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
308 if (fQLongLCMSNotCalc)
312 Double_t beta = fPzSum/fESum;
313 fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
314 fQLongLCMSNotCalc = kFALSE;
318 /************************************************************************/
320 Double_t AliAODPair::GetQtLCMS()
322 //returns Q transverse CMS longitudionally co-moving
325 fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
326 fQtLCMSNotCalc = kFALSE;
330 /************************************************************************/
332 Double_t AliAODPair::GetQt()
334 //returns Q transverse CMS longitudionally co-moving
340 Double_t dotprod = fPxSum*fPxDiff + fPySum*fPyDiff + fPzSum*fPzDiff;
341 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
342 klen = TMath::Sqrt(klen);
343 Double_t qlen = fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff;
344 qlen = TMath::Sqrt(qlen);
345 Double_t tmp = klen*qlen;
352 Double_t cosopenangle = dotprod/tmp;
353 Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
355 fQt = sinopenangle*qlen;
360 /************************************************************************/
362 Double_t AliAODPair::GetKt()
364 //calculates the evarage momentum of the pair
368 fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
373 /************************************************************************/
375 Double_t AliAODPair::GetKStar()
377 //calculates invariant velocity difference
382 Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
383 Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
386 Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
387 Error("GetKStar","Particle1:");
389 Error("GetKStar","Particle2:");
391 Error("GetKStar","");
394 fKStarNotCalc = kFALSE;
397 Double_t pinv = TMath::Sqrt(mtrans - ptrans);
399 Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
406 Info("GetKStar","Sqrt of negative number q = %f",q);
407 Error("GetKStar","Particle1:");
409 Error("GetKStar","Particle2:");
412 fKStarNotCalc = kFALSE;
418 fKStarNotCalc = kFALSE;
422 /************************************************************************/
423 Double_t AliAODPair::GetKStarOut()
425 CalculateKStarComp();
428 /************************************************************************/
429 Double_t AliAODPair::GetKStarSide()
431 CalculateKStarComp();
434 /************************************************************************/
435 Double_t AliAODPair::GetKStarLong()
437 CalculateKStarComp();
440 /************************************************************************/
442 Double_t AliAODPair::GetQInv()
445 //warning for non-id particles you want to use 2*KStar
449 fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
450 fQInvNotCalc = kFALSE;
454 /************************************************************************/
456 Double_t AliAODPair::GetGammaToLCMS()
458 //calculates gamma factor of the boost to LCMS
459 if(fGammaLCMSNotCalc)
462 Double_t beta = fPzSum/fESum;
463 fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
464 fGammaLCMSNotCalc = kFALSE;
468 /************************************************************************/
470 Double_t AliAODPair::GetGammaToTransverse()
472 //calculates gamma factor of the boost to LCMS
473 Double_t beta = 2.0*GetKt() / GetMt();
474 Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
478 /************************************************************************/
480 Double_t AliAODPair::GetMt()
482 //Calculates transverse mass of the pair
486 fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
491 /************************************************************************/
493 Double_t AliAODPair::GetAvarageDistance()
495 //returns and buffers avarage distance between two tracks calculated
496 // out of track points (see AliAODTrackPoints class)
498 if (fAvarageDistanceNotCalc)
500 fAvarageDistance = AvDistance();
501 fAvarageDistanceNotCalc = kFALSE;
503 return fAvarageDistance;
505 /************************************************************************/
507 Double_t AliAODPair::AvDistance()
509 //returns avarage distance between two tracks in range
510 //as defined in Track-Points of AliVAODParticle
511 //returns negative value if error uccured f.g. tracks do not have track-points
512 AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
514 {//it could be simulated pair
515 // Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
519 AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
522 // Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
526 return tpts1->AvarageDistance(*tpts2);
528 /************************************************************************/
530 Double_t AliAODPair::GetR()
532 //Returns distance between particles vertexes in thir CMS
536 Double_t vxDiff = fPart1->Vx() - fPart2->Vx();
537 Double_t vyDiff = fPart1->Vy() - fPart2->Vy();
538 Double_t vzDiff = fPart1->Vz() - fPart2->Vz();
540 Double_t r = TMath::Sqrt( vxDiff*vxDiff + vyDiff*vyDiff + vzDiff*vzDiff );
544 /************************************************************************/
546 Double_t AliAODPair::GetRStar()
548 //Returns distance between particles vertexes in thir CMS
553 Double_t klen = fPxSum*fPxSum + fPySum*fPySum + fPzSum*fPzSum;
554 klen = TMath::Sqrt(klen);
556 Double_t aBeta = klen/fESum;
557 Double_t aGamma = 1.0/TMath::Sqrt(1.0 - aBeta*aBeta);
560 Double_t alpha = -TMath::ATan2(fPySum,fPzSum);
561 Double_t beta = TMath::ATan2(fPxSum,TMath::Hypot(fPySum,fPzSum));
563 Double_t sinalpha = TMath::Sin(alpha);
564 Double_t cosalpha = TMath::Cos(alpha);
565 Double_t sinbeta = TMath::Sin(beta);
566 Double_t cosbeta = TMath::Cos(beta);
568 Double_t v1xP = fPart1->Vx();
569 Double_t v2xP = fPart2->Vx();
570 Double_t v1yP = fPart1->Vy()*cosalpha + fPart1->Vz()*sinalpha;
571 Double_t v2yP = fPart2->Vy()*cosalpha + fPart2->Vz()*sinalpha;
572 Double_t v1zP =-fPart1->Vy()*sinalpha + fPart1->Vz()*cosalpha;
573 Double_t v2zP =-fPart2->Vy()*sinalpha + fPart2->Vz()*cosalpha;
576 ///////////////////////////////////////////////////
578 // Double_t p1yP = fPart1->Py()*cosalpha + fPart1->Pz()*sinalpha;
579 // Double_t p2yP = fPart2->Py()*cosalpha + fPart2->Pz()*sinalpha;
581 // Double_t p1zP =-fPart1->Py()*sinalpha + fPart1->Pz()*cosalpha;
582 // Double_t p2zP =-fPart2->Py()*sinalpha + fPart2->Pz()*cosalpha;
585 // Double_t p1x = fPart1->Px()*cosbeta - p1zP*sinbeta;
586 // Double_t p2x = fPart2->Px()*cosbeta - p2zP*sinbeta;
587 // Double_t p1z = fPart1->Px()*sinbeta + p1zP*cosbeta;
588 // Double_t p2z = fPart2->Px()*sinbeta + p2zP*cosbeta;
590 // Info("","%f %f %f",p1yP,p2yP,p1yP+p2yP);
591 // Info("","%f %f %f",p1x,p2x,p1x+p2x);
593 // Info("","%f %f ",p1x+p2x,p1yP+p2yP);
595 ///////////////////////////////////////////////////
598 Double_t v1x = v1xP*cosbeta - v1zP*sinbeta;
599 Double_t v2x = v2xP*cosbeta - v2zP*sinbeta;
602 Double_t v1z = v1xP*sinbeta + v1zP*cosbeta;
603 Double_t v2z = v2xP*sinbeta + v2zP*cosbeta;
606 Double_t v1zB=aGamma*(v1z-aBeta*fPart1->T());
607 Double_t v2zB=aGamma*(v2z-aBeta*fPart2->T());
611 Double_t dx = v1x - v2x;
612 Double_t dy = v1y - v2y;
613 Double_t dz = v1zB - v2zB;
615 Double_t rstar = TMath::Sqrt( dx*dx + dy*dy + dz*dz);
619 /************************************************************************/
621 void AliAODPair::MirrorSecond()
623 //makes local copy of the second particle and mirrors their momenta
624 //for its deletion is responsible who calls this method
625 fPart2 = (AliVAODParticle*)fPart2->Clone();
626 fPart2->SetMomentum(-fPart2->Px(),-fPart2->Py(),-fPart2->Pz(),fPart2->E());
629 /************************************************************************/
631 void AliAODPair::DeleteSecond()
633 //Deletes second particle
638 void AliAODPair::Print()
640 if (fPart1) fPart1->Print();
641 if (fPart2) fPart2->Print();
643 Info("Print","GetKStar() %f",GetKStar());
644 Info("Print","GetKt() %f",GetKt() );
645 Info("Print","QInv %f", GetQInv() );
646 Info("Print","GetQOutLCMS() %f",GetQOutLCMS() );
647 Info("Print","GetQSideLCMS %f",GetQSideLCMS() );
648 Info("Print","GetQLongLCMS() %f", GetQLongLCMS());
649 Info("Print","GetDeltaTheta() %f", GetDeltaTheta());
650 Info("Print","GetDeltaPhi() %f", GetDeltaPhi());