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