Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackPair.cxx
CommitLineData
820b4d9e 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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//====================================================================================================================================================
17//
18// Description of an ALICE muon forward track pair, i.e. a pair of AliMuonForwardTrack objects
19//
20// Contact author: antonio.uras@cern.ch
21//
22//====================================================================================================================================================
23
24#include "AliLog.h"
25#include "AliMUONTrackParam.h"
26#include "TParticle.h"
27#include "AliMuonForwardTrack.h"
28#include "AliMUONTrackExtrap.h"
29#include "TClonesArray.h"
30#include "TDatabasePDG.h"
31#include "TLorentzVector.h"
32#include "TRandom.h"
33#include "AliMuonForwardTrackPair.h"
34
35ClassImp(AliMuonForwardTrackPair)
36
37//====================================================================================================================================================
38
39AliMuonForwardTrackPair::AliMuonForwardTrackPair():
40 TObject(),
d4643a10 41 fMuonForwardTracks(0),
bcaf50eb 42 fKinemMC(0,0,0,0),
43 fKinem(0,0,0,0),
a2b7dc2a 44 fIsKinemSet(kFALSE),
45 fXPointOfClosestApproach(9999),
46 fYPointOfClosestApproach(9999),
47 fZPointOfClosestApproach(9999)
820b4d9e 48{
49
50 // default constructor
51
52 fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
274c2dce 53 fMuonForwardTracks -> SetOwner(kTRUE);
54
820b4d9e 55}
56
57//====================================================================================================================================================
58
59AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack *track0, AliMuonForwardTrack *track1):
60 TObject(),
d4643a10 61 fMuonForwardTracks(0),
bcaf50eb 62 fKinemMC(0,0,0,0),
63 fKinem(0,0,0,0),
a2b7dc2a 64 fIsKinemSet(kFALSE),
65 fXPointOfClosestApproach(9999),
66 fYPointOfClosestApproach(9999),
67 fZPointOfClosestApproach(9999)
820b4d9e 68{
69
70 fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
274c2dce 71 fMuonForwardTracks->SetOwner(kTRUE);
820b4d9e 72 new ((*fMuonForwardTracks)[0]) AliMuonForwardTrack(*track0);
73 new ((*fMuonForwardTracks)[1]) AliMuonForwardTrack(*track1);
74
d4643a10 75 SetKinemMC();
a2b7dc2a 76 SetPointOfClosestApproach();
d4643a10 77
820b4d9e 78}
79
80//====================================================================================================================================================
81
82AliMuonForwardTrackPair::AliMuonForwardTrackPair(const AliMuonForwardTrackPair& trackPair):
83 TObject(trackPair),
274c2dce 84 fMuonForwardTracks(NULL),
bcaf50eb 85 fKinemMC(trackPair.fKinemMC),
86 fKinem(trackPair.fKinem),
a2b7dc2a 87 fIsKinemSet(trackPair.fIsKinemSet),
88 fXPointOfClosestApproach(trackPair.fXPointOfClosestApproach),
89 fYPointOfClosestApproach(trackPair.fYPointOfClosestApproach),
90 fZPointOfClosestApproach(trackPair.fZPointOfClosestApproach)
820b4d9e 91{
92
93 // copy constructor
274c2dce 94 fMuonForwardTracks = new TClonesArray(*(trackPair.fMuonForwardTracks));
95 fMuonForwardTracks -> SetOwner(kTRUE);
96
820b4d9e 97}
98
99//====================================================================================================================================================
100
101AliMuonForwardTrackPair& AliMuonForwardTrackPair::operator=(const AliMuonForwardTrackPair& trackPair) {
102
103 // Asignment operator
104
105 // check assignement to self
106 if (this == &trackPair) return *this;
107
108 // base class assignement
109 AliMuonForwardTrackPair::operator=(trackPair);
110
111 // clear memory
274c2dce 112 Clear("");
113
114 fMuonForwardTracks = new TClonesArray(*(trackPair.fMuonForwardTracks));
115 fMuonForwardTracks->SetOwner(kTRUE);
820b4d9e 116
bcaf50eb 117 fKinemMC = trackPair.fKinemMC;
118 fKinem = trackPair.fKinem;
119 fIsKinemSet = trackPair.fIsKinemSet;
a2b7dc2a 120 fXPointOfClosestApproach = trackPair.fXPointOfClosestApproach;
121 fYPointOfClosestApproach = trackPair.fYPointOfClosestApproach;
122 fZPointOfClosestApproach = trackPair.fZPointOfClosestApproach;
820b4d9e 123
124 return *this;
125
126}
127
128//====================================================================================================================================================
129
820b4d9e 130Double_t AliMuonForwardTrackPair::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
131
132 Double_t weightedOffset[2]={0};
133
a6f7b1e5 134 weightedOffset[0] = GetTrack(0)->GetWeightedOffset(x, y, z);
135 weightedOffset[1] = GetTrack(1)->GetWeightedOffset(x, y, z);
820b4d9e 136
137 Double_t weightedOffsetDimuon = TMath::Sqrt(0.5 * (weightedOffset[0]*weightedOffset[0] + weightedOffset[1]*weightedOffset[1]));
138
139 return weightedOffsetDimuon;
140
141}
142
cabfe25b 143//====================================================================================================================================================
144
145Double_t AliMuonForwardTrackPair::GetWeightedOffsetAtPCA() {
146
147 return GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
148
149}
150
151//====================================================================================================================================================
152
153Double_t AliMuonForwardTrackPair::GetPCAQuality() {
154
155 Double_t wOffset_1 = GetTrack(0)->GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
156 Double_t wOffset_2 = GetTrack(1)->GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
157
158 Double_t f1 = TMath::Exp(-0.5 * wOffset_1);
159 Double_t f2 = TMath::Exp(-0.5 * wOffset_2);
160
161 if (f1+f2 > 0.) return (2.*f1*f2)/(f1+f2);
162 else return 0.;
163
164}
165
820b4d9e 166//====================================================================================================================================================
167
820b4d9e 168Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Double_t z, Int_t nClusters) {
169
170 Int_t idCluster[2] = {0};
171 if (nClusters>0) {
a6f7b1e5 172 idCluster[0] = GetTrack(0)->GetNMUONClusters() - nClusters;
173 idCluster[1] = GetTrack(1)->GetNMUONClusters() - nClusters;
820b4d9e 174 }
175 if (idCluster[0]<0) idCluster[0] = 0;
176 if (idCluster[1]<0) idCluster[1] = 0;
177
a6f7b1e5 178 AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMUONCluster(idCluster[0]);
179 AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMUONCluster(idCluster[1]);
820b4d9e 180
d4643a10 181 AliDebug(2, Form("MUON before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5 182 param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
d4643a10 183
184 AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
820b4d9e 185 AliMUONTrackExtrap::ExtrapToVertex(param0, x, y, z, 0., 0.); // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...)
d4643a10 186 AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
820b4d9e 187 AliMUONTrackExtrap::ExtrapToVertex(param1, x, y, z, 0., 0.); // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...)
188
d4643a10 189 AliDebug(2, Form("MUON after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5 190 param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
d4643a10 191
820b4d9e 192 Double_t momentum[2] = {0};
193
194 momentum[0] = param0->P();
195 momentum[1] = param1->P();
196
197 Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
198
199 TLorentzVector dimu;
200
201 dimu.SetE(TMath::Sqrt(mMu*mMu + momentum[0]*momentum[0]) + TMath::Sqrt(mMu*mMu + momentum[1]*momentum[1]));
202
203 dimu.SetPx(param0->Px() + param1->Px());
204 dimu.SetPy(param0->Py() + param1->Py());
205 dimu.SetPz(param0->Pz() + param1->Pz());
206
207 return dimu.M();
208
209}
210
211//====================================================================================================================================================
212
d4643a10 213void AliMuonForwardTrackPair::SetKinemMC() {
7230691d 214
a6f7b1e5 215 if ( !(GetTrack(0)->GetMCTrackRef()) || !(GetTrack(1)->GetMCTrackRef()) ) return;
bcaf50eb 216
d4643a10 217 AliDebug(2, Form("MC: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5 218 GetTrack(0)->GetMCTrackRef()->Px(), GetTrack(0)->GetMCTrackRef()->Py(), GetTrack(0)->GetMCTrackRef()->Pz(),
219 GetTrack(1)->GetMCTrackRef()->Px(), GetTrack(1)->GetMCTrackRef()->Py(), GetTrack(1)->GetMCTrackRef()->Pz()));
220
221 fKinemMC.SetE(GetTrack(0)->GetMCTrackRef()->Energy() + GetTrack(1)->GetMCTrackRef()->Energy());
d4643a10 222
a6f7b1e5 223 fKinemMC.SetPx(GetTrack(0)->GetMCTrackRef()->Px() + GetTrack(1)->GetMCTrackRef()->Px());
224 fKinemMC.SetPy(GetTrack(0)->GetMCTrackRef()->Py() + GetTrack(1)->GetMCTrackRef()->Py());
225 fKinemMC.SetPz(GetTrack(0)->GetMCTrackRef()->Pz() + GetTrack(1)->GetMCTrackRef()->Pz());
820b4d9e 226
227}
228
229//====================================================================================================================================================
bcaf50eb 230
231void AliMuonForwardTrackPair::SetKinem(Double_t z, Int_t nClusters) {
232
233// if (!fMuonForwardTracks) return kFALSE;
234// if (!fMuonForwardTracks->At(0) || !fMuonForwardTracks->At(1)) return kFALSE;
235
236 Int_t idCluster[2] = {0};
237 if (nClusters>0) {
a6f7b1e5 238 idCluster[0] = GetTrack(0)->GetNMFTClusters() - nClusters;
239 idCluster[1] = GetTrack(1)->GetNMFTClusters() - nClusters;
bcaf50eb 240 }
241 if (idCluster[0]<0) idCluster[0] = 0;
242 if (idCluster[1]<0) idCluster[1] = 0;
243
244 Double_t momentum[2] = {0};
245
a6f7b1e5 246 AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(idCluster[0]);
247 AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(idCluster[1]);
bcaf50eb 248
249 AliDebug(2, Form("MFT before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5 250 param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
bcaf50eb 251
252 if (TMath::Abs(z)<1e6) {
253 AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
254 AliMUONTrackExtrap::ExtrapToZCov(param0, z);
255 AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
256 AliMUONTrackExtrap::ExtrapToZCov(param1, z);
257 }
258
259 AliDebug(2, Form("MFT after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5 260 param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
bcaf50eb 261
262 momentum[0] = (param0->P());
263 momentum[1] = (param1->P());
264
265 Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
266
267 fKinem.SetE(TMath::Sqrt(mMu*mMu + momentum[0]*momentum[0]) + TMath::Sqrt(mMu*mMu + momentum[1]*momentum[1]));
268 fKinem.SetPx(param0->Px() + param1->Px());
269 fKinem.SetPy(param0->Py() + param1->Py());
270 fKinem.SetPz(param0->Pz() + param1->Pz());
271
272 fIsKinemSet = kTRUE;
273
a2b7dc2a 274}
275
276//====================================================================================================================================================
277
278void AliMuonForwardTrackPair::SetPointOfClosestApproach() {
279
a6f7b1e5 280 AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(0);
281 AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(0);
a2b7dc2a 282
283 Double_t step = 1.; // in cm
284 Double_t startPoint = 0.;
bcaf50eb 285
a2b7dc2a 286 Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
287
288 for (Int_t i=0; i<3; i++) {
289 AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
290 AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
291 Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
292 Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
293 r[i] = TMath::Sqrt(dX*dX + dY*dY);
294 }
295
296 Double_t researchDirection=0.;
297
298 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.; // towards z positive
299 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.; // towards z negative
300 else if (r[0]<r[1] && r[1]>r[2]) {
301 AliError("Point of closest approach cannot be found for dimuon (no minima)");
302 return;
303 }
304
305 while (TMath::Abs(researchDirection)>0.5) {
306 if (researchDirection>0.) {
307 z[0] = z[1];
308 z[1] = z[2];
309 z[2] = z[1]+researchDirection*step;
310 }
311 else {
312 z[2] = z[1];
313 z[1] = z[0];
314 z[0] = z[1]+researchDirection*step;
315 }
316 if (TMath::Abs(z[0])>900.) {
317 AliError("Point of closest approach cannot be found for dimuon (no minima)");
318 return;
319 }
320 for (Int_t i=0; i<3; i++) {
321 AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
322 AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
323 Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
324 Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
325 r[i] = TMath::Sqrt(dX*dX + dY*dY);
326 }
327 researchDirection=0.;
328 if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.; // towards z positive
329 else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.; // towards z negative
330 }
331
332 AliDebug(2,"Minimum region has been found");
333
334 step *= 0.5;
335 while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
336 z[0] = z[1]-step;
337 z[2] = z[1]+step;
338 for (Int_t i=0; i<3; i++) {
339 AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
340 AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
341 Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
342 Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
343 r[i] = TMath::Sqrt(dX*dX + dY*dY);
344 }
345 if (r[0]<r[1]) z[1] = z[0];
346 else if (r[2]<r[1]) z[1] = z[2];
347 else step *= 0.5;
348 }
349
350 fZPointOfClosestApproach = z[1];
351 AliMUONTrackExtrap::ExtrapToZCov(param0, fZPointOfClosestApproach);
352 AliMUONTrackExtrap::ExtrapToZCov(param1, fZPointOfClosestApproach);
353 fXPointOfClosestApproach = 0.5*(param0->GetNonBendingCoor() + param1->GetNonBendingCoor());
354 fYPointOfClosestApproach = 0.5*(param0->GetBendingCoor() + param1->GetBendingCoor());
355
356 AliDebug(2,Form("Point of Closest Approach: (%f, %f, %f)",fXPointOfClosestApproach,fYPointOfClosestApproach,fZPointOfClosestApproach));
357
bcaf50eb 358}
359
360//====================================================================================================================================================
820b4d9e 361
d4643a10 362Bool_t AliMuonForwardTrackPair::IsResonance() {
363
364 Bool_t result = kFALSE;
820b4d9e 365
d4643a10 366 Int_t labelMC[2] = {0};
367 Int_t codePDG[2] = {0};
368
369 for (Int_t iTrack=0; iTrack<2; iTrack++) {
a6f7b1e5 370 labelMC[iTrack] = GetTrack(iTrack)->GetParentMCLabel(0);
371 codePDG[iTrack] = GetTrack(iTrack)->GetParentPDGCode(0);
d4643a10 372 }
820b4d9e 373
d4643a10 374 AliDebug(1, Form("Muons' mothers: (%d, %d)", labelMC[0], labelMC[1]));
820b4d9e 375
d4643a10 376 if (labelMC[0]==labelMC[1] && codePDG[0]==codePDG[1] && (codePDG[0]== 113 ||
7cb5b5b9 377 codePDG[0]== 221 ||
d4643a10 378 codePDG[0]== 223 ||
7cb5b5b9 379 codePDG[0]== 331 ||
d4643a10 380 codePDG[0]== 333 ||
381 codePDG[0]== 443 ||
382 codePDG[0]==100443 ||
383 codePDG[0]== 553 ||
384 codePDG[0]==100553 ) ) result = kTRUE;
820b4d9e 385
d4643a10 386 if (result) AliDebug(1, Form("Pair is a resonance %d", codePDG[0]));
820b4d9e 387
d4643a10 388 return result;
820b4d9e 389
d4643a10 390}
820b4d9e 391
d4643a10 392//====================================================================================================================================================
45c2d1ca 393