Bug in storage manager making possible to delete permanent event fixed
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackPair.cxx
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
35 ClassImp(AliMuonForwardTrackPair)
36
37 //====================================================================================================================================================
38
39 AliMuonForwardTrackPair::AliMuonForwardTrackPair():
40   TObject(),
41   fMuonForwardTracks(0),
42   fKinemMC(0,0,0,0),
43   fKinem(0,0,0,0),
44   fIsKinemSet(kFALSE),
45   fXPointOfClosestApproach(9999),
46   fYPointOfClosestApproach(9999),
47   fZPointOfClosestApproach(9999)
48 {
49
50   // default constructor
51
52   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
53   fMuonForwardTracks -> SetOwner(kTRUE);
54   
55 }
56
57 //====================================================================================================================================================
58
59 AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack *track0, AliMuonForwardTrack *track1):
60   TObject(),
61   fMuonForwardTracks(0),
62   fKinemMC(0,0,0,0),
63   fKinem(0,0,0,0),
64   fIsKinemSet(kFALSE),
65   fXPointOfClosestApproach(9999),
66   fYPointOfClosestApproach(9999),
67   fZPointOfClosestApproach(9999)
68 {
69
70   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
71   fMuonForwardTracks->SetOwner(kTRUE);
72   new ((*fMuonForwardTracks)[0]) AliMuonForwardTrack(*track0);
73   new ((*fMuonForwardTracks)[1]) AliMuonForwardTrack(*track1);
74
75   SetKinemMC();
76   SetPointOfClosestApproach();
77
78 }
79
80 //====================================================================================================================================================
81
82 AliMuonForwardTrackPair::AliMuonForwardTrackPair(const AliMuonForwardTrackPair& trackPair): 
83   TObject(trackPair),
84   fMuonForwardTracks(NULL),
85   fKinemMC(trackPair.fKinemMC),
86   fKinem(trackPair.fKinem),
87   fIsKinemSet(trackPair.fIsKinemSet),
88   fXPointOfClosestApproach(trackPair.fXPointOfClosestApproach),
89   fYPointOfClosestApproach(trackPair.fYPointOfClosestApproach),
90   fZPointOfClosestApproach(trackPair.fZPointOfClosestApproach)
91 {
92
93   // copy constructor
94   fMuonForwardTracks = new TClonesArray(*(trackPair.fMuonForwardTracks));
95   fMuonForwardTracks -> SetOwner(kTRUE);
96
97 }
98
99 //====================================================================================================================================================
100
101 AliMuonForwardTrackPair& 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
112   Clear("");
113   
114   fMuonForwardTracks      = new TClonesArray(*(trackPair.fMuonForwardTracks));
115   fMuonForwardTracks->SetOwner(kTRUE);
116   
117   fKinemMC = trackPair.fKinemMC;
118   fKinem = trackPair.fKinem;
119   fIsKinemSet = trackPair.fIsKinemSet;
120   fXPointOfClosestApproach = trackPair.fXPointOfClosestApproach;
121   fYPointOfClosestApproach = trackPair.fYPointOfClosestApproach;
122   fZPointOfClosestApproach = trackPair.fZPointOfClosestApproach;
123
124   return *this;
125
126 }
127
128 //====================================================================================================================================================
129
130 Double_t AliMuonForwardTrackPair::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
131
132   Double_t weightedOffset[2]={0};
133
134   weightedOffset[0] = GetTrack(0)->GetWeightedOffset(x, y, z);
135   weightedOffset[1] = GetTrack(1)->GetWeightedOffset(x, y, z);
136
137   Double_t weightedOffsetDimuon = TMath::Sqrt(0.5 * (weightedOffset[0]*weightedOffset[0] + weightedOffset[1]*weightedOffset[1]));
138
139   return weightedOffsetDimuon;
140
141 }
142
143 //====================================================================================================================================================
144
145 Double_t AliMuonForwardTrackPair::GetWeightedOffsetAtPCA() {
146
147   return GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
148
149 }
150
151 //====================================================================================================================================================
152
153 Double_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
166 //====================================================================================================================================================
167
168 Double_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) {
172     idCluster[0] = GetTrack(0)->GetNMUONClusters() - nClusters;
173     idCluster[1] = GetTrack(1)->GetNMUONClusters() - nClusters;
174   }
175   if (idCluster[0]<0) idCluster[0] = 0;
176   if (idCluster[1]<0) idCluster[1] = 0;
177
178   AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMUONCluster(idCluster[0]);
179   AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMUONCluster(idCluster[1]);
180
181   AliDebug(2, Form("MUON before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
182                    param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
183
184   AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
185   AliMUONTrackExtrap::ExtrapToVertex(param0, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
186   AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
187   AliMUONTrackExtrap::ExtrapToVertex(param1, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
188
189   AliDebug(2, Form("MUON after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
190                    param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
191
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
213 void AliMuonForwardTrackPair::SetKinemMC() {
214
215   if ( !(GetTrack(0)->GetMCTrackRef()) || !(GetTrack(1)->GetMCTrackRef()) ) return;
216
217   AliDebug(2, Form("MC: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
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());
222   
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());
226
227 }
228
229 //====================================================================================================================================================
230
231 void 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) {
238     idCluster[0] = GetTrack(0)->GetNMFTClusters() - nClusters;
239     idCluster[1] = GetTrack(1)->GetNMFTClusters() - nClusters;
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   
246   AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(idCluster[0]);
247   AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(idCluster[1]);
248
249   AliDebug(2, Form("MFT before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
250                    param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
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)", 
260                    param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
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
274 }
275
276 //====================================================================================================================================================
277
278 void AliMuonForwardTrackPair::SetPointOfClosestApproach() {
279   
280   AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(0);
281   AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(0);
282   
283   Double_t step = 1.;  // in cm
284   Double_t startPoint = 0.;
285
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   
358 }
359
360 //====================================================================================================================================================
361
362 Bool_t AliMuonForwardTrackPair::IsResonance() {
363
364   Bool_t result = kFALSE;
365
366   Int_t labelMC[2] = {0};
367   Int_t codePDG[2] = {0};
368   
369   for (Int_t iTrack=0; iTrack<2; iTrack++) {
370     labelMC[iTrack] = GetTrack(iTrack)->GetParentMCLabel(0);
371     codePDG[iTrack] = GetTrack(iTrack)->GetParentPDGCode(0);
372   }
373
374   AliDebug(1, Form("Muons' mothers: (%d, %d)", labelMC[0], labelMC[1]));
375
376   if (labelMC[0]==labelMC[1] && codePDG[0]==codePDG[1] && (codePDG[0]==   113 ||
377                                                            codePDG[0]==   221 ||
378                                                            codePDG[0]==   223 ||
379                                                            codePDG[0]==   331 ||
380                                                            codePDG[0]==   333 ||
381                                                            codePDG[0]==   443 ||
382                                                            codePDG[0]==100443 ||
383                                                            codePDG[0]==   553 ||
384                                                            codePDG[0]==100553 ) ) result = kTRUE;
385
386   if (result) AliDebug(1, Form("Pair is a resonance %d", codePDG[0]));
387
388   return result; 
389
390 }
391
392 //====================================================================================================================================================
393