]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMuonForwardTrack.cxx
Temporary bug fixed in AliMFTTrackerMU
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrack.cxx
index 0a5897f015d7dc26301754ecba2221a6f70030e8..46851a10c2bfb3b80f7923f5668930d1a778c4ab 100644 (file)
 #include "TMatrixD.h"
 #include "TParticle.h"
 #include "AliMuonForwardTrack.h"
+#include "AliMFTConstants.h"
+#include "TLorentzVector.h"
+#include "TDatabasePDG.h"
+#include "AliMUONConstants.h"
 
 ClassImp(AliMuonForwardTrack)
 
@@ -40,12 +44,23 @@ AliMuonForwardTrack::AliMuonForwardTrack():
   AliMUONTrack(),
   fMUONTrack(0),
   fMCTrackRef(0),
-  fMFTClusters(0)
+  fMFTClusters(0),
+  fNWrongClustersMC(-1),
+  fTrackMCId(-1),
+  fKinem(0,0,0,0),
+  fParamCovMatrix(5,5),
+  fRAtAbsorberEnd(0)
 {
 
   // default constructor
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = -1;
+    fParentPDGCode[iParent] =  0;
+  }
   fMFTClusters = new TClonesArray("AliMFTCluster");
+  fMFTClusters -> SetOwner(kTRUE);
 
 }
 
@@ -55,12 +70,22 @@ AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
   AliMUONTrack(),
   fMUONTrack(0),
   fMCTrackRef(0),
-  fMFTClusters(0)
+  fMFTClusters(0),
+  fNWrongClustersMC(-1),
+  fTrackMCId(-1),
+  fKinem(0,0,0,0),
+  fParamCovMatrix(5,5),
+  fRAtAbsorberEnd(0)
 {
 
   SetMUONTrack(MUONTrack);
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = -1;
+    fParentPDGCode[iParent] =  0;
+  }
   fMFTClusters = new TClonesArray("AliMFTCluster");
+  fMFTClusters -> SetOwner(kTRUE);
 
 }
 
@@ -68,13 +93,26 @@ AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
 
 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
   AliMUONTrack(track),
-  fMUONTrack(track.fMUONTrack),
-  fMCTrackRef(track.fMCTrackRef),
-  fMFTClusters(track.fMFTClusters)
+  fMUONTrack(0x0),
+  fMCTrackRef(0x0),
+  fMFTClusters(0x0),
+  fNWrongClustersMC(track.fNWrongClustersMC),
+  fTrackMCId(track.fTrackMCId),
+  fKinem(track.fKinem),
+  fParamCovMatrix(track.fParamCovMatrix),
+  fRAtAbsorberEnd(track.fRAtAbsorberEnd)
 {
 
   // copy constructor
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
+  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
+  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
+  fMFTClusters->SetOwner(kTRUE);
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
+    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
+  }
   
 }
 
@@ -91,13 +129,23 @@ AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& t
   AliMUONTrack::operator=(track);
   
   // clear memory
-  Clear();
+  Clear("");
   
-  fMUONTrack    = track.fMUONTrack;
-  fMCTrackRef   = track.fMCTrackRef;
-  fMFTClusters  = track.fMFTClusters;
-
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
+  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
+  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
+  fMFTClusters->SetOwner(kTRUE);
+  fNWrongClustersMC = track.fNWrongClustersMC;
+  fTrackMCId        = track.fTrackMCId;
+  fKinem            = track.fKinem;
+  fParamCovMatrix   = track.fParamCovMatrix;
+  fRAtAbsorberEnd   = track.fRAtAbsorberEnd;
+
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
+    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
+  }
   
   return *this;
 
@@ -105,6 +153,29 @@ AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& t
 
 //====================================================================================================================================================
 
+void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
+
+  // Clear arrays
+  fMFTClusters -> Delete(); 
+  delete fMFTClusters; fMFTClusters = 0x0;
+  delete fMUONTrack;   fMUONTrack   = 0x0;
+  delete fMCTrackRef;  fMCTrackRef  = 0x0;
+  
+}
+
+//====================================================================================================================================================
+
+AliMuonForwardTrack::~AliMuonForwardTrack() {
+
+  delete fMUONTrack;
+  delete fMCTrackRef;
+  fMFTClusters -> Delete();
+  delete fMFTClusters;
+  
+}
+
+//====================================================================================================================================================
+
 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
 
   if (fMUONTrack) {
@@ -113,6 +184,7 @@ void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
   }
 
   fMUONTrack = MUONTrack;
+  SetGlobalChi2(fMUONTrack->GetGlobalChi2());
   
 }
 
@@ -133,8 +205,8 @@ void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
 
 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
 
-  AliDebug(1, Form("fMFTClusters=%p has %d entries", fMFTClusters, fMFTClusters->GetEntries()));
-  Int_t iMFTCluster = fMFTClusters->GetEntries();
+  AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
+  Int_t iMFTCluster = this->fMFTClusters->GetEntries();
   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
           mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
@@ -142,15 +214,18 @@ void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackPara
   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
-  AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries",GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries()));
   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
+  AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
   mftCluster.SetLocalChi2(chi2Kalman);
   mftCluster.SetTrackChi2(newGlobalChi2);
-  new ((*fMFTClusters)[iMFTCluster]) AliMFTCluster(mftCluster);
+  new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
+  AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
+                  GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
   SetGlobalChi2(newGlobalChi2);
+  AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
@@ -210,12 +285,12 @@ AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
 
   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
-    AliError("Invalid MFT cluster index. NULL pointer will be returned");
+    AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
     return NULL;
   }
 
   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
-  AliMFTCluster *mftCluster = (AliMFTCluster*) fMFTClusters->At(trackParam->GetUniqueID());
+  AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
 
   return mftCluster;
 
@@ -361,3 +436,306 @@ Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
 
 //====================================================================================================================================================
 
+Double_t AliMuonForwardTrack::GetDCA(Double_t x, Double_t y, Double_t z) {
+
+  // Distance of Closest Approach, according to the standard MUON terminology. Actually, the offset of the track w.r.t. the primary vertex,
+  // where the extrapolation of the track DOES NOT include the MFT information
+
+  AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
+
+  AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&param, z);
+  Double_t dX = param.GetNonBendingCoor() - x;
+  Double_t dY = param.GetBendingCoor()    - y;
+  return TMath::Sqrt(dX*dX + dY*dY);
+
+}
+
+//====================================================================================================================================================
+
+Double_t AliMuonForwardTrack::GetMomentumSpectrometer(Double_t z) {
+
+  // Momentum of the track at the primary vertex plane, where the extrapolation of the track DOES NOT include the MFT information
+
+  AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
+
+  AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&param, z);
+  return param.P();
+
+}
+
+//====================================================================================================================================================
+
+Double_t AliMuonForwardTrack::GetThetaAbs() {
+
+  // it is the angle defined by the imaginary line goingo from the vertex to the exit point of the track at the end of the hadron absorber
+  
+  Double_t z = AliMUONConstants::AbsZEnd();
+
+  AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
+
+  AliMUONTrackExtrap::ExtrapToZ(&param, z);
+  Double_t x = param.GetNonBendingCoor();
+  Double_t y = param.GetBendingCoor();
+
+  return 180. +TMath::ATan(TMath::Sqrt(x*x + y*y)/z)*TMath::RadToDeg();
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromDirectResonance() {
+
+  if (!IsMuon()) return kFALSE;
+
+  for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGResonance(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-resonance particle
+
+  return kTRUE; 
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromChainResonance() {
+
+  if (!IsMuon()) return kFALSE;
+
+  if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
+
+  if (!IsPDGResonance(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a resonance
+
+  for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGResonance(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-resonance particle
+
+  return kFALSE;
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromDirectCharm() {
+
+  if (!IsMuon()) return kFALSE;
+
+  for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGCharm(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-charmed particle
+
+  return kTRUE; 
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromChainCharm() {
+
+  if (!IsMuon()) return kFALSE;
+
+  if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
+
+  if (!IsPDGCharm(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a charmed hadron
+
+  for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGCharm(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-charmed particle
+
+  return kFALSE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromDirectBeauty() {
+
+  if (!IsMuon()) return kFALSE;
+
+  for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGBeauty(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-beauty particle
+
+  return kTRUE; 
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFromChainBeauty() {
+
+  if (!IsMuon()) return kFALSE;
+
+  if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
+
+  if (!IsPDGBeauty(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a beauty hadron
+
+  for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGBeauty(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-beauty particle
+
+  return kFALSE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsMuon() {
+
+  if (IsFake()) return kFALSE;
+
+  if (TMath::Abs(fMCTrackRef->GetPdgCode())==13) return kTRUE;
+
+  return kFALSE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsFake() {
+
+  if (!fMCTrackRef) return kTRUE;
+
+  return kFALSE;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsPDGResonance(Int_t pdg) {
+
+  // if (pdg<10) return kFALSE; 
+  // Int_t id = pdg%100000; 
+  // return (!((id-id%10)%110));
+
+  Bool_t result = kFALSE;
+
+  if ( pdg ==    113 ||
+       pdg ==    221 ||
+       pdg ==    223 ||
+       pdg ==    331 ||
+       pdg ==    333 ||
+       pdg ==    443 ||
+       pdg == 100443 ||
+       pdg ==    553 ||
+       pdg == 100553 ) result = kTRUE;
+  
+  return result; 
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsPDGCharm(Int_t pdg) {
+
+  Bool_t result = kFALSE;
+
+  if ( TMath::Abs(pdg) ==   411 ||
+       TMath::Abs(pdg) ==   421 ||
+       TMath::Abs(pdg) == 10411 ||
+       TMath::Abs(pdg) == 10421 ||
+       TMath::Abs(pdg) ==   413 ||
+       TMath::Abs(pdg) ==   423 ||
+       TMath::Abs(pdg) == 10413 ||
+       TMath::Abs(pdg) == 10423 ||
+       TMath::Abs(pdg) == 20413 ||
+       TMath::Abs(pdg) == 20423 ||
+       TMath::Abs(pdg) ==   415 ||
+       TMath::Abs(pdg) ==   425 ||
+       TMath::Abs(pdg) ==   431 ||
+       TMath::Abs(pdg) == 10431 ||
+       TMath::Abs(pdg) ==   433 ||
+       TMath::Abs(pdg) == 10433 ||
+       TMath::Abs(pdg) == 20433 ||
+       TMath::Abs(pdg) ==   435 ) result = kTRUE;
+  
+  return result; 
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsPDGBeauty(Int_t pdg) {
+
+  Bool_t result = kFALSE;
+
+  if ( TMath::Abs(pdg) ==   511 ||
+       TMath::Abs(pdg) ==   521 ||
+       TMath::Abs(pdg) == 10511 ||
+       TMath::Abs(pdg) == 10521 ||
+       TMath::Abs(pdg) ==   513 ||
+       TMath::Abs(pdg) ==   523 ||
+       TMath::Abs(pdg) == 10513 ||
+       TMath::Abs(pdg) == 10523 ||
+       TMath::Abs(pdg) == 20513 ||
+       TMath::Abs(pdg) == 20523 ||
+       TMath::Abs(pdg) ==   515 ||
+       TMath::Abs(pdg) ==   525 ||
+       TMath::Abs(pdg) ==   531 ||
+       TMath::Abs(pdg) == 10531 ||
+       TMath::Abs(pdg) ==   533 ||
+       TMath::Abs(pdg) == 10533 ||
+       TMath::Abs(pdg) == 20533 ||
+       TMath::Abs(pdg) ==   535 ||
+       TMath::Abs(pdg) ==   541 ||
+       TMath::Abs(pdg) == 10541 ||
+       TMath::Abs(pdg) ==   543 ||
+       TMath::Abs(pdg) == 10543 ||
+       TMath::Abs(pdg) == 20543 ||
+       TMath::Abs(pdg) ==   545 ) result = kTRUE;
+  
+  return result; 
+  
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMuonForwardTrack::IsMuonFromBackground() {
+
+  Bool_t result = kFALSE;
+
+  if (!IsMuon()) return result;
+
+  if (!IsFromDirectResonance() && !IsFromChainResonance() && !IsFromDirectCharm() && !IsFromChainCharm() && !IsFromDirectBeauty() && !IsFromChainBeauty()) result = kTRUE;
+
+  if (result) AliDebug(1, Form("Muon comes from a background source %d", GetParentPDGCode(0)));
+
+  return result;
+
+}
+
+//====================================================================================================================================================
+
+void AliMuonForwardTrack::EvalKinem(Double_t z) {
+
+  AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
+  AliMUONTrackExtrap::ExtrapToZCov(param, z);
+
+  Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
+  Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
+  fKinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
+
+  TMatrixD cov(5,5);
+  fParamCovMatrix = param->GetCovariances();
+
+}
+
+//====================================================================================================================================================
+
+Int_t AliMuonForwardTrack::GetFirstMotherID() {
+
+  Int_t motherLevel = 0; 
+  Int_t motherMCLabel = GetParentMCLabel(motherLevel); 
+  while (motherMCLabel >= 0) { 
+    motherLevel++; 
+    motherMCLabel = GetParentMCLabel(motherLevel); 
+  }
+  return motherLevel-1;
+
+}
+
+//====================================================================================================================================================
+
+void AliMuonForwardTrack::PrintHistory() {
+
+  if (IsFake()) printf("Track is a fake MUON\n");
+  else {
+    TString history = "";
+    for (Int_t i=GetFirstMotherID(); i>=0; i--) {
+      history += TDatabasePDG::Instance()->GetParticle(GetParentPDGCode(i))->GetName();
+      history += " -> ";
+    }
+    history += TDatabasePDG::Instance()->GetParticle(fMCTrackRef->GetPdgCode())->GetName();
+    printf("%s\n",history.Data());
+  }
+
+}
+
+//====================================================================================================================================================