[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackPair.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//====================================================================================================================================================
//
//      Description of an ALICE muon forward track pair, i.e. a pair of AliMuonForwardTrack objects
//
//      Contact author: antonio.uras@cern.ch
//
//====================================================================================================================================================

#include "AliLog.h"
#include "AliMUONTrackParam.h"
#include "TParticle.h"
#include "AliMuonForwardTrack.h"
#include "AliMUONTrackExtrap.h"
#include "TClonesArray.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"
#include "TRandom.h"
#include "AliMuonForwardTrackPair.h"

ClassImp(AliMuonForwardTrackPair)

//====================================================================================================================================================

AliMuonForwardTrackPair::AliMuonForwardTrackPair():
  TObject(),
  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),
  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(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);

}

//====================================================================================================================================================

AliMuonForwardTrackPair& AliMuonForwardTrackPair::operator=(const AliMuonForwardTrackPair& trackPair) {

  // Asignment operator

  // check assignement to self
  if (this == &trackPair) return *this;

  // base class assignement
  AliMuonForwardTrackPair::operator=(trackPair);
  
  // clear memory
  Clear("");
  
  fMuonForwardTracks      = new TClonesArray(*(trackPair.fMuonForwardTracks));
  fMuonForwardTracks->SetOwner(kTRUE);
  
  fKinemMC = trackPair.fKinemMC;
  fKinem = trackPair.fKinem;
  fIsKinemSet = trackPair.fIsKinemSet;
  fXPointOfClosestApproach = trackPair.fXPointOfClosestApproach;
  fYPointOfClosestApproach = trackPair.fYPointOfClosestApproach;
  fZPointOfClosestApproach = trackPair.fZPointOfClosestApproach;

  return *this;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrackPair::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {

  Double_t weightedOffset[2]={0};

  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]));

  return weightedOffsetDimuon;

}

//====================================================================================================================================================

Double_t AliMuonForwardTrackPair::GetWeightedOffsetAtPCA() {

  return GetWeightedOffset(fXPointOfClosestApproach, fYPointOfClosestApproach, fZPointOfClosestApproach);

}

//====================================================================================================================================================

Double_t AliMuonForwardTrackPair::GetPCAQuality() {

  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.;
  
}

//====================================================================================================================================================

Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Double_t z, Int_t nClusters) {

  Int_t idCluster[2] = {0};
  if (nClusters>0) {
    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 = 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();
  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(param0->Px() + param1->Px());
  dimu.SetPy(param0->Py() + param1->Py());
  dimu.SetPz(param0->Pz() + param1->Pz());

  return dimu.M();

}

//====================================================================================================================================================

void AliMuonForwardTrackPair::SetKinemMC() {

  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()));

  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());

}

//====================================================================================================================================================

void AliMuonForwardTrackPair::SetKinem(Double_t z, Int_t nClusters) {

//   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; 

}

//====================================================================================================================================================
Commit	Line	Data
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
	35	ClassImp(AliMuonForwardTrackPair)
	36
	37	//====================================================================================================================================================
	38
	39	AliMuonForwardTrackPair::AliMuonForwardTrackPair():
	40	TObject(),
d4643a10	41	fMuonForwardTracks(0),
bcaf50eb	42	fKinemMC(0,0,0,0),
bcaf50eb	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);
274c2dce	54
820b4d9e	55	}
	56
	57	//====================================================================================================================================================
	58
	59	AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack track0, AliMuonForwardTrack track1):
	60	TObject(),
d4643a10	61	fMuonForwardTracks(0),
bcaf50eb	62	fKinemMC(0,0,0,0),
bcaf50eb	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
	82	AliMuonForwardTrackPair::AliMuonForwardTrackPair(const AliMuonForwardTrackPair& trackPair):
	83	TObject(trackPair),
274c2dce	84	fMuonForwardTracks(NULL),
bcaf50eb	85	fKinemMC(trackPair.fKinemMC),
bcaf50eb	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
	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
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	130	Double_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);
a6f7b1e5	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
	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.f1f2)/(f1+f2);
	162	else return 0.;
	163
	164	}
	165
820b4d9e	166	//====================================================================================================================================================
820b4d9e	167
820b4d9e	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) {
a6f7b1e5	172	idCluster[0] = GetTrack(0)->GetNMUONClusters() - nClusters;
a6f7b1e5	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]);
a6f7b1e5	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
d4643a10	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(...)
820b4d9e	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(mMumMu + momentum[0]momentum[0]) + TMath::Sqrt(mMumMu + 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	213	void 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
	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) {
a6f7b1e5	238	idCluster[0] = GetTrack(0)->GetNMFTClusters() - nClusters;
a6f7b1e5	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]);
a6f7b1e5	247	AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(idCluster[1]);
bcaf50eb	248
bcaf50eb	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(mMumMu + momentum[0]momentum[0]) + TMath::Sqrt(mMumMu + 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
	278	void AliMuonForwardTrackPair::SetPointOfClosestApproach() {
	279
a6f7b1e5	280	AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(0);
a6f7b1e5	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(dXdX + dYdY);
	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(dXdX + dYdY);
	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(dXdX + dYdY);
	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	362	Bool_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);
a6f7b1e5	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	//====================================================================================================================================================