]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMuonForwardTrackPair.cxx
Update to FindFASTJET.cmake; now accepts also for -DFASTJET= option
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackPair.cxx
index ac23f1c6f76bb9edc4684e21903d6ad55ff9a811..4bc7684d4b237ededf9de2761359cc534034c302 100644 (file)
@@ -38,38 +38,62 @@ ClassImp(AliMuonForwardTrackPair)
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair():
   TObject(),
-  fMuonForwardTracks(0)
+  fMuonForwardTracks(0),
+  fKinemMC(0,0,0,0),
+  fKinem(0,0,0,0),
+  fIsKinemSet(kFALSE),
+  fXPointOfClosestApproach(9999),
+  fYPointOfClosestApproach(9999),
+  fZPointOfClosestApproach(9999)
 {
 
   // default constructor
 
   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
-
+  fMuonForwardTracks -> SetOwner(kTRUE);
+  
 }
 
 //====================================================================================================================================================
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack *track0, AliMuonForwardTrack *track1):
   TObject(),
-  fMuonForwardTracks(0)
+  fMuonForwardTracks(0),
+  fKinemMC(0,0,0,0),
+  fKinem(0,0,0,0),
+  fIsKinemSet(kFALSE),
+  fXPointOfClosestApproach(9999),
+  fYPointOfClosestApproach(9999),
+  fZPointOfClosestApproach(9999)
 {
 
   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
-
+  fMuonForwardTracks->SetOwner(kTRUE);
   new ((*fMuonForwardTracks)[0]) AliMuonForwardTrack(*track0);
   new ((*fMuonForwardTracks)[1]) AliMuonForwardTrack(*track1);
 
+  SetKinemMC();
+  SetPointOfClosestApproach();
+
 }
 
 //====================================================================================================================================================
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair(const AliMuonForwardTrackPair& trackPair): 
   TObject(trackPair),
-  fMuonForwardTracks(trackPair.fMuonForwardTracks)
+  fMuonForwardTracks(NULL),
+  fKinemMC(trackPair.fKinemMC),
+  fKinem(trackPair.fKinem),
+  fIsKinemSet(trackPair.fIsKinemSet),
+  fXPointOfClosestApproach(trackPair.fXPointOfClosestApproach),
+  fYPointOfClosestApproach(trackPair.fYPointOfClosestApproach),
+  fZPointOfClosestApproach(trackPair.fZPointOfClosestApproach)
 {
 
   // copy constructor
-  
+  fMuonForwardTracks = new TClonesArray(*(trackPair.fMuonForwardTracks));
+  fMuonForwardTracks -> SetOwner(kTRUE);
+
 }
 
 //====================================================================================================================================================
@@ -85,9 +109,17 @@ AliMuonForwardTrackPair& AliMuonForwardTrackPair::operator=(const AliMuonForward
   AliMuonForwardTrackPair::operator=(trackPair);
   
   // clear memory
-  Clear();
+  Clear("");
+  
+  fMuonForwardTracks      = new TClonesArray(*(trackPair.fMuonForwardTracks));
+  fMuonForwardTracks->SetOwner(kTRUE);
   
-  fMuonForwardTracks = trackPair.fMuonForwardTracks;
+  fKinemMC = trackPair.fKinemMC;
+  fKinem = trackPair.fKinem;
+  fIsKinemSet = trackPair.fIsKinemSet;
+  fXPointOfClosestApproach = trackPair.fXPointOfClosestApproach;
+  fYPointOfClosestApproach = trackPair.fYPointOfClosestApproach;
+  fZPointOfClosestApproach = trackPair.fZPointOfClosestApproach;
 
   return *this;
 
@@ -95,22 +127,12 @@ AliMuonForwardTrackPair& AliMuonForwardTrackPair::operator=(const AliMuonForward
 
 //====================================================================================================================================================
 
-void AliMuonForwardTrackPair::SetTrack(Int_t iTrack, AliMuonForwardTrack *track) {
-
-  if (iTrack==0 || iTrack==1) {
-    new ((*fMuonForwardTracks)[iTrack]) AliMuonForwardTrack(*track);
-  }
-
-}
-
-//====================================================================================================================================================
-
 Double_t AliMuonForwardTrackPair::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
 
   Double_t weightedOffset[2]={0};
 
-  weightedOffset[0] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetWeightedOffset(x, y, z);
-  weightedOffset[1] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetWeightedOffset(x, y, z);
+  weightedOffset[0] = GetTrack(0)->GetWeightedOffset(x, y, z);
+  weightedOffset[1] = GetTrack(1)->GetWeightedOffset(x, y, z);
 
   Double_t weightedOffsetDimuon = TMath::Sqrt(0.5 * (weightedOffset[0]*weightedOffset[0] + weightedOffset[1]*weightedOffset[1]));
 
@@ -120,46 +142,25 @@ Double_t AliMuonForwardTrackPair::GetWeightedOffset(Double_t x, Double_t y, Doub
 
 //====================================================================================================================================================
 
-Double_t AliMuonForwardTrackPair::GetMass(Double_t z, Int_t nClusters) {
-
-  Int_t idCluster[2] = {0};
-  if (nClusters>0) {
-    idCluster[0] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetNMFTClusters() - nClusters;
-    idCluster[1] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetNMFTClusters() - nClusters;
-  }
-  if (idCluster[0]<0) idCluster[0] = 0;
-  if (idCluster[1]<0) idCluster[1] = 0;
+Double_t AliMuonForwardTrackPair::GetWeightedOffsetAtPCA() {
 
-  Double_t momentum[2] = {0};
-  
-  AliMUONTrackParam *param0 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0]);
-  AliMUONTrackParam *param1 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1]);
+  return GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
 
-  if (TMath::Abs(z)<1e6) {
-    AliMUONTrackExtrap::ExtrapToZCov(param0, z);
-    AliMUONTrackExtrap::ExtrapToZCov(param1, z);
-  }
-
-  momentum[0] = (param0->P());
-  momentum[1] = (param1->P());
-
-  Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
-
-  TLorentzVector dimu;
-
-  dimu.SetE(TMath::Sqrt(mMu*mMu + momentum[0]*momentum[0]) + TMath::Sqrt(mMu*mMu + momentum[1]*momentum[1]));
-
-  dimu.SetPx(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Px()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Px());
+}
 
-  dimu.SetPy(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Py()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Py());
+//====================================================================================================================================================
 
-  dimu.SetPz(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Pz()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Pz());
+Double_t AliMuonForwardTrackPair::GetPCAQuality() {
 
-  return dimu.M();
+  Double_t wOffset_1 = GetTrack(0)->GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
+  Double_t wOffset_2 = GetTrack(1)->GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);
+  
+  Double_t f1 = TMath::Exp(-0.5 * wOffset_1);
+  Double_t f2 = TMath::Exp(-0.5 * wOffset_2);
 
+  if (f1+f2 > 0.) return (2.*f1*f2)/(f1+f2); 
+  else return 0.;
+  
 }
 
 //====================================================================================================================================================
@@ -168,18 +169,26 @@ Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Doub
 
   Int_t idCluster[2] = {0};
   if (nClusters>0) {
-    idCluster[0] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetNMUONClusters() - nClusters;
-    idCluster[1] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetNMUONClusters() - nClusters;
+    idCluster[0] = GetTrack(0)->GetNMUONClusters() - nClusters;
+    idCluster[1] = GetTrack(1)->GetNMUONClusters() - nClusters;
   }
   if (idCluster[0]<0) idCluster[0] = 0;
   if (idCluster[1]<0) idCluster[1] = 0;
 
-  AliMUONTrackParam *param0 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMUONCluster(idCluster[0]);
-  AliMUONTrackParam *param1 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMUONCluster(idCluster[1]);
+  AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMUONCluster(idCluster[0]);
+  AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMUONCluster(idCluster[1]);
+
+  AliDebug(2, Form("MUON before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
 
+  AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
   AliMUONTrackExtrap::ExtrapToVertex(param0, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
+  AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
   AliMUONTrackExtrap::ExtrapToVertex(param1, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
 
+  AliDebug(2, Form("MUON after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
+
   Double_t momentum[2] = {0}; 
 
   momentum[0] = param0->P();
@@ -201,32 +210,184 @@ Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Doub
 
 //====================================================================================================================================================
 
-Double_t AliMuonForwardTrackPair::GetMassMC() {
+void AliMuonForwardTrackPair::SetKinemMC() {
 
-  TLorentzVector dimu;
+  if ( !(GetTrack(0)->GetMCTrackRef()) || !(GetTrack(1)->GetMCTrackRef()) ) return;
+
+  AliDebug(2, Form("MC: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  GetTrack(0)->GetMCTrackRef()->Px(), GetTrack(0)->GetMCTrackRef()->Py(), GetTrack(0)->GetMCTrackRef()->Pz(),
+                  GetTrack(1)->GetMCTrackRef()->Px(), GetTrack(1)->GetMCTrackRef()->Py(), GetTrack(1)->GetMCTrackRef()->Pz()));
 
-  dimu.SetE(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Energy() +
-           ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Energy());
+  fKinemMC.SetE(GetTrack(0)->GetMCTrackRef()->Energy() + GetTrack(1)->GetMCTrackRef()->Energy());
+  
+  fKinemMC.SetPx(GetTrack(0)->GetMCTrackRef()->Px() + GetTrack(1)->GetMCTrackRef()->Px());
+  fKinemMC.SetPy(GetTrack(0)->GetMCTrackRef()->Py() + GetTrack(1)->GetMCTrackRef()->Py());
+  fKinemMC.SetPz(GetTrack(0)->GetMCTrackRef()->Pz() + GetTrack(1)->GetMCTrackRef()->Pz());
 
-  dimu.SetPx(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Px() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Px());
+}
 
-  dimu.SetPy(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Py() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Py());
+//====================================================================================================================================================
 
-  dimu.SetPz(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Pz() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Pz());
+void AliMuonForwardTrackPair::SetKinem(Double_t z, Int_t nClusters) {
 
-  return dimu.M();
+//   if (!fMuonForwardTracks) return kFALSE;
+//   if (!fMuonForwardTracks->At(0) || !fMuonForwardTracks->At(1)) return kFALSE;
+
+  Int_t idCluster[2] = {0};
+  if (nClusters>0) {
+    idCluster[0] = GetTrack(0)->GetNMFTClusters() - nClusters;
+    idCluster[1] = GetTrack(1)->GetNMFTClusters() - nClusters;
+  }
+  if (idCluster[0]<0) idCluster[0] = 0;
+  if (idCluster[1]<0) idCluster[1] = 0;
+
+  Double_t momentum[2] = {0};
+  
+  AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(idCluster[0]);
+  AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(idCluster[1]);
+
+  AliDebug(2, Form("MFT before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
+
+  if (TMath::Abs(z)<1e6) {
+    AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
+    AliMUONTrackExtrap::ExtrapToZCov(param0, z);
+    AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
+    AliMUONTrackExtrap::ExtrapToZCov(param1, z);
+  }
+
+  AliDebug(2, Form("MFT after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
+
+  momentum[0] = (param0->P());
+  momentum[1] = (param1->P());
+
+  Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
+
+  fKinem.SetE(TMath::Sqrt(mMu*mMu + momentum[0]*momentum[0]) + TMath::Sqrt(mMu*mMu + momentum[1]*momentum[1]));
+  fKinem.SetPx(param0->Px() + param1->Px());
+  fKinem.SetPy(param0->Py() + param1->Py());
+  fKinem.SetPz(param0->Pz() + param1->Pz());
+
+  fIsKinemSet = kTRUE;
 
 }
 
 //====================================================================================================================================================
 
+void AliMuonForwardTrackPair::SetPointOfClosestApproach() {
+  
+  AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(0);
+  AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(0);
+  
+  Double_t step = 1.;  // in cm
+  Double_t startPoint = 0.;
 
+  Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
+  
+  for (Int_t i=0; i<3; i++) {
+    AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
+    AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
+    Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+    Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+    r[i] = TMath::Sqrt(dX*dX + dY*dY);
+  }
+  
+  Double_t researchDirection=0.;
+  
+  if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.;   // towards z positive
+  else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.;   // towards z negative
+  else if (r[0]<r[1] && r[1]>r[2]) { 
+    AliError("Point of closest approach cannot be found for dimuon (no minima)");
+    return;
+  }
+  
+  while (TMath::Abs(researchDirection)>0.5) {
+    if (researchDirection>0.) {
+      z[0] = z[1];
+      z[1] = z[2];
+      z[2] = z[1]+researchDirection*step;
+    }
+    else {
+      z[2] = z[1];
+      z[1] = z[0];
+      z[0] = z[1]+researchDirection*step;
+    }
+    if (TMath::Abs(z[0])>900.) {
+      AliError("Point of closest approach cannot be found for dimuon (no minima)");
+      return;
+    }
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
+      AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
+      Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+      Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+      r[i] = TMath::Sqrt(dX*dX + dY*dY);
+    }
+    researchDirection=0.;
+    if      (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.;   // towards z positive
+    else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.;   // towards z negative
+  }
+  
+  AliDebug(2,"Minimum region has been found");
+  
+  step *= 0.5;
+  while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
+    z[0] = z[1]-step;
+    z[2] = z[1]+step;
+    for (Int_t i=0; i<3; i++) {
+      AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
+      AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
+      Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
+      Double_t dY = param0->GetBendingCoor()    - param1->GetBendingCoor();
+      r[i] = TMath::Sqrt(dX*dX + dY*dY);
+    }
+    if      (r[0]<r[1]) z[1] = z[0];
+    else if (r[2]<r[1]) z[1] = z[2];
+    else step *= 0.5;
+  }
+  
+  fZPointOfClosestApproach = z[1];
+  AliMUONTrackExtrap::ExtrapToZCov(param0, fZPointOfClosestApproach);
+  AliMUONTrackExtrap::ExtrapToZCov(param1, fZPointOfClosestApproach);  
+  fXPointOfClosestApproach = 0.5*(param0->GetNonBendingCoor() + param1->GetNonBendingCoor());
+  fYPointOfClosestApproach = 0.5*(param0->GetBendingCoor()    + param1->GetBendingCoor());
+  
+  AliDebug(2,Form("Point of Closest Approach: (%f, %f, %f)",fXPointOfClosestApproach,fYPointOfClosestApproach,fZPointOfClosestApproach));
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrackPair::IsResonance() {
+
+  Bool_t result = kFALSE;
 
+  Int_t labelMC[2] = {0};
+  Int_t codePDG[2] = {0};
+  
+  for (Int_t iTrack=0; iTrack<2; iTrack++) {
+    labelMC[iTrack] = GetTrack(iTrack)->GetParentMCLabel(0);
+    codePDG[iTrack] = GetTrack(iTrack)->GetParentPDGCode(0);
+  }
+
+  AliDebug(1, Form("Muons' mothers: (%d, %d)", labelMC[0], labelMC[1]));
 
+  if (labelMC[0]==labelMC[1] && codePDG[0]==codePDG[1] && (codePDG[0]==   113 ||
+                                                          codePDG[0]==   221 ||
+                                                          codePDG[0]==   223 ||
+                                                          codePDG[0]==   331 ||
+                                                          codePDG[0]==   333 ||
+                                                          codePDG[0]==   443 ||
+                                                          codePDG[0]==100443 ||
+                                                          codePDG[0]==   553 ||
+                                                          codePDG[0]==100553 ) ) result = kTRUE;
 
+  if (result) AliDebug(1, Form("Pair is a resonance %d", codePDG[0]));
 
+  return result; 
 
+}
+
+//====================================================================================================================================================