[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::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
	143	//====================================================================================================================================================
	144
820b4d9e	145	Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Double_t z, Int_t nClusters) {
	146
	147	Int_t idCluster[2] = {0};
	148	if (nClusters>0) {
a6f7b1e5	149	idCluster[0] = GetTrack(0)->GetNMUONClusters() - nClusters;
a6f7b1e5	150	idCluster[1] = GetTrack(1)->GetNMUONClusters() - nClusters;
820b4d9e	151	}
	152	if (idCluster[0]<0) idCluster[0] = 0;
	153	if (idCluster[1]<0) idCluster[1] = 0;
	154
a6f7b1e5	155	AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMUONCluster(idCluster[0]);
a6f7b1e5	156	AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMUONCluster(idCluster[1]);
820b4d9e	157
d4643a10	158	AliDebug(2, Form("MUON before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5	159	param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
d4643a10	160
d4643a10	161	AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
820b4d9e	162	AliMUONTrackExtrap::ExtrapToVertex(param0, x, y, z, 0., 0.); // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...)
d4643a10	163	AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
820b4d9e	164	AliMUONTrackExtrap::ExtrapToVertex(param1, x, y, z, 0., 0.); // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...)
820b4d9e	165
d4643a10	166	AliDebug(2, Form("MUON after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5	167	param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
d4643a10	168
820b4d9e	169	Double_t momentum[2] = {0};
	170
	171	momentum[0] = param0->P();
	172	momentum[1] = param1->P();
	173
	174	Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
	175
	176	TLorentzVector dimu;
	177
	178	dimu.SetE(TMath::Sqrt(mMumMu + momentum[0]momentum[0]) + TMath::Sqrt(mMumMu + momentum[1]momentum[1]));
	179
	180	dimu.SetPx(param0->Px() + param1->Px());
	181	dimu.SetPy(param0->Py() + param1->Py());
	182	dimu.SetPz(param0->Pz() + param1->Pz());
	183
	184	return dimu.M();
	185
	186	}
	187
	188	//====================================================================================================================================================
	189
d4643a10	190	void AliMuonForwardTrackPair::SetKinemMC() {
7230691d	191
a6f7b1e5	192	if ( !(GetTrack(0)->GetMCTrackRef()) \|\| !(GetTrack(1)->GetMCTrackRef()) ) return;
bcaf50eb	193
d4643a10	194	AliDebug(2, Form("MC: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5	195	GetTrack(0)->GetMCTrackRef()->Px(), GetTrack(0)->GetMCTrackRef()->Py(), GetTrack(0)->GetMCTrackRef()->Pz(),
	196	GetTrack(1)->GetMCTrackRef()->Px(), GetTrack(1)->GetMCTrackRef()->Py(), GetTrack(1)->GetMCTrackRef()->Pz()));
	197
	198	fKinemMC.SetE(GetTrack(0)->GetMCTrackRef()->Energy() + GetTrack(1)->GetMCTrackRef()->Energy());
d4643a10	199
a6f7b1e5	200	fKinemMC.SetPx(GetTrack(0)->GetMCTrackRef()->Px() + GetTrack(1)->GetMCTrackRef()->Px());
	201	fKinemMC.SetPy(GetTrack(0)->GetMCTrackRef()->Py() + GetTrack(1)->GetMCTrackRef()->Py());
	202	fKinemMC.SetPz(GetTrack(0)->GetMCTrackRef()->Pz() + GetTrack(1)->GetMCTrackRef()->Pz());
820b4d9e	203
	204	}
	205
	206	//====================================================================================================================================================
bcaf50eb	207
	208	void AliMuonForwardTrackPair::SetKinem(Double_t z, Int_t nClusters) {
	209
	210	// if (!fMuonForwardTracks) return kFALSE;
	211	// if (!fMuonForwardTracks->At(0) \|\| !fMuonForwardTracks->At(1)) return kFALSE;
	212
	213	Int_t idCluster[2] = {0};
	214	if (nClusters>0) {
a6f7b1e5	215	idCluster[0] = GetTrack(0)->GetNMFTClusters() - nClusters;
a6f7b1e5	216	idCluster[1] = GetTrack(1)->GetNMFTClusters() - nClusters;
bcaf50eb	217	}
	218	if (idCluster[0]<0) idCluster[0] = 0;
	219	if (idCluster[1]<0) idCluster[1] = 0;
	220
	221	Double_t momentum[2] = {0};
	222
a6f7b1e5	223	AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(idCluster[0]);
a6f7b1e5	224	AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(idCluster[1]);
bcaf50eb	225
bcaf50eb	226	AliDebug(2, Form("MFT before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5	227	param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
bcaf50eb	228
	229	if (TMath::Abs(z)<1e6) {
	230	AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
	231	AliMUONTrackExtrap::ExtrapToZCov(param0, z);
	232	AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
	233	AliMUONTrackExtrap::ExtrapToZCov(param1, z);
	234	}
	235
	236	AliDebug(2, Form("MFT after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)",
a6f7b1e5	237	param0->Px(), param0->Py(), param0->Pz(), param1->Px(), param1->Py(), param1->Pz()));
bcaf50eb	238
	239	momentum[0] = (param0->P());
	240	momentum[1] = (param1->P());
	241
	242	Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
	243
	244	fKinem.SetE(TMath::Sqrt(mMumMu + momentum[0]momentum[0]) + TMath::Sqrt(mMumMu + momentum[1]momentum[1]));
	245	fKinem.SetPx(param0->Px() + param1->Px());
	246	fKinem.SetPy(param0->Py() + param1->Py());
	247	fKinem.SetPz(param0->Pz() + param1->Pz());
	248
	249	fIsKinemSet = kTRUE;
	250
a2b7dc2a	251	}
	252
	253	//====================================================================================================================================================
	254
	255	void AliMuonForwardTrackPair::SetPointOfClosestApproach() {
	256
a6f7b1e5	257	AliMUONTrackParam *param0 = GetTrack(0)->GetTrackParamAtMFTCluster(0);
a6f7b1e5	258	AliMUONTrackParam *param1 = GetTrack(1)->GetTrackParamAtMFTCluster(0);
a2b7dc2a	259
	260	Double_t step = 1.; // in cm
	261	Double_t startPoint = 0.;
bcaf50eb	262
a2b7dc2a	263	Double_t r[3]={0}, z[3]={startPoint, startPoint+step, startPoint+2*step};
	264
	265	for (Int_t i=0; i<3; i++) {
	266	AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
	267	AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
	268	Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
	269	Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
	270	r[i] = TMath::Sqrt(dXdX + dYdY);
	271	}
	272
	273	Double_t researchDirection=0.;
	274
	275	if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.; // towards z positive
	276	else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.; // towards z negative
	277	else if (r[0]<r[1] && r[1]>r[2]) {
	278	AliError("Point of closest approach cannot be found for dimuon (no minima)");
	279	return;
	280	}
	281
	282	while (TMath::Abs(researchDirection)>0.5) {
	283	if (researchDirection>0.) {
	284	z[0] = z[1];
	285	z[1] = z[2];
	286	z[2] = z[1]+researchDirection*step;
	287	}
	288	else {
	289	z[2] = z[1];
	290	z[1] = z[0];
	291	z[0] = z[1]+researchDirection*step;
	292	}
	293	if (TMath::Abs(z[0])>900.) {
	294	AliError("Point of closest approach cannot be found for dimuon (no minima)");
	295	return;
	296	}
	297	for (Int_t i=0; i<3; i++) {
	298	AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
	299	AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
	300	Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
	301	Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
	302	r[i] = TMath::Sqrt(dXdX + dYdY);
	303	}
	304	researchDirection=0.;
	305	if (r[0]>r[1] && r[1]>r[2]) researchDirection = +1.; // towards z positive
	306	else if (r[0]<r[1] && r[1]<r[2]) researchDirection = -1.; // towards z negative
	307	}
	308
	309	AliDebug(2,"Minimum region has been found");
	310
	311	step *= 0.5;
	312	while (step>AliMFTConstants::fPrecisionPointOfClosestApproach) {
	313	z[0] = z[1]-step;
	314	z[2] = z[1]+step;
	315	for (Int_t i=0; i<3; i++) {
	316	AliMUONTrackExtrap::ExtrapToZCov(param0, z[i]);
	317	AliMUONTrackExtrap::ExtrapToZCov(param1, z[i]);
	318	Double_t dX = param0->GetNonBendingCoor() - param1->GetNonBendingCoor();
	319	Double_t dY = param0->GetBendingCoor() - param1->GetBendingCoor();
	320	r[i] = TMath::Sqrt(dXdX + dYdY);
	321	}
	322	if (r[0]<r[1]) z[1] = z[0];
	323	else if (r[2]<r[1]) z[1] = z[2];
	324	else step *= 0.5;
	325	}
	326
327	fZPointOfClosestApproach = z[1];
328	AliMUONTrackExtrap::ExtrapToZCov(param0, fZPointOfClosestApproach);
329	AliMUONTrackExtrap::ExtrapToZCov(param1, fZPointOfClosestApproach);
330	fXPointOfClosestApproach = 0.5*(param0->GetNonBendingCoor() + param1->GetNonBendingCoor());
331	fYPointOfClosestApproach = 0.5*(param0->GetBendingCoor() + param1->GetBendingCoor());
332
333	AliDebug(2,Form("Point of Closest Approach: (%f, %f, %f)",fXPointOfClosestApproach,fYPointOfClosestApproach,fZPointOfClosestApproach));
334
bcaf50eb	335	}
	336
	337	//====================================================================================================================================================
820b4d9e	338
d4643a10	339	Bool_t AliMuonForwardTrackPair::IsResonance() {
	340
	341	Bool_t result = kFALSE;
820b4d9e	342
d4643a10	343	Int_t labelMC[2] = {0};
	344	Int_t codePDG[2] = {0};
	345
	346	for (Int_t iTrack=0; iTrack<2; iTrack++) {
a6f7b1e5	347	labelMC[iTrack] = GetTrack(iTrack)->GetParentMCLabel(0);
a6f7b1e5	348	codePDG[iTrack] = GetTrack(iTrack)->GetParentPDGCode(0);
d4643a10	349	}
820b4d9e	350
d4643a10	351	AliDebug(1, Form("Muons' mothers: (%d, %d)", labelMC[0], labelMC[1]));
820b4d9e	352
d4643a10	353	if (labelMC[0]==labelMC[1] && codePDG[0]==codePDG[1] && (codePDG[0]== 113 \|\|
7cb5b5b9	354	codePDG[0]== 221 \|\|
d4643a10	355	codePDG[0]== 223 \|\|
7cb5b5b9	356	codePDG[0]== 331 \|\|
d4643a10	357	codePDG[0]== 333 \|\|
	358	codePDG[0]== 443 \|\|
	359	codePDG[0]==100443 \|\|
	360	codePDG[0]== 553 \|\|
	361	codePDG[0]==100553 ) ) result = kTRUE;
820b4d9e	362
d4643a10	363	if (result) AliDebug(1, Form("Pair is a resonance %d", codePDG[0]));
820b4d9e	364
d4643a10	365	return result;
820b4d9e	366
d4643a10	367	}
820b4d9e	368
d4643a10	369	//====================================================================================================================================================