[u/mrichter/AliRoot.git] / MUON / AliMUONTrackLight.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *      SigmaEffect_thetadegrees                                          *
 * 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 purpeateose. It is      *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-----------------------------------------------------------------------------
// Compact information for the muon generated tracks in the MUON arm 
// useful at the last stage of the analysis chain
// provides a link between the reconstructed track and the generated particle 
// stores kinematical information at gen. and rec. level and 
// the decay history of the muon, allowing the identification of the 
// mother process 
// 
// To be used together with AliMUONPairLight
//
// This class was prepared by INFN Cagliari, July 2006
// (authors: H.Woehri, A.de Falco)
//-----------------------------------------------------------------------------

// 13 Nov 2007:
// Added a temporary fix to FindRefTrack to be able to handle reconstructed tracks
// generated from ESD muon track information. The problem is that the ESD data at
// the moment only contains the first hit on chamber 1. Hopefully in the near future
// this will be fixed and all hit information will be available.
//  - Artur Szostak <artursz@iafrica.com>

#include "AliMUONTrackLight.h"
#include "AliMUONTrack.h"
#include "AliMUONConstants.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrackParam.h"

#include "AliESDMuonTrack.h"
#include "AliStack.h"
#include "AliLog.h"

#include "TDatabasePDG.h"
#include "TParticle.h"
#include "TString.h"

#include <cstdio>

ClassImp(AliMUONTrackLight) 

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

AliMUONTrackLight::AliMUONTrackLight() 
  : TObject(), 
    fPrec(), 
    fIsTriggered(kFALSE),
    fCharge(-999), 
    fChi2(-1), 
    fCentr(-1),
    fPgen(), 
    fTrackPythiaLine(-999),
    fTrackPDGCode(-999),
    fOscillation(kFALSE), 
    fNParents(0),
    fWeight(1)    
{
  /// default constructor
  fPgen.SetPxPyPzE(0.,0.,0.,0.); 
  fPrec.SetPxPyPzE(0.,0.,0.,0.); 
  for (Int_t i=0; i<3; i++) fXYZ[i]=-999; 
  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
    fParentPDGCode[npar] = -1; 
    fParentPythiaLine[npar] = -1;
  }
  for (Int_t i = 0; i < 4; i++){
    fQuarkPDGCode[i] = -1; 
    fQuarkPythiaLine[i] = -1; 
  }
}

//============================================
AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy) 
  : TObject(muonCopy), 
    fPrec(muonCopy.fPrec), 
    fIsTriggered(muonCopy.fIsTriggered),
    fCharge(muonCopy.fCharge), 
    fChi2(muonCopy.fChi2), 
    fCentr(muonCopy.fCentr),
    fPgen(muonCopy.fPgen), 
    fTrackPythiaLine(muonCopy.fTrackPythiaLine),
    fTrackPDGCode(muonCopy.fTrackPDGCode),
    fOscillation(muonCopy.fOscillation), 
    fNParents(muonCopy.fNParents),
    fWeight(muonCopy.fWeight)
{
  /// copy constructor
  for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i]; 
  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
    fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar]; 
    fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
  }
  for (Int_t i = 0; i < 4; i++){
    fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i]; 
    fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i]; 
  }
}

//============================================
AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
  : TObject(), 
    fPrec(), 
    fIsTriggered(kFALSE),
    fCharge(-999), 
    fChi2(-1),
    fCentr(-1),
    fPgen(), 
    fTrackPythiaLine(-999),
    fTrackPDGCode(-999),
    fOscillation(kFALSE), 
    fNParents(0),
    fWeight(1)
{ 
  /// constructor
  fPgen.SetPxPyPzE(0.,0.,0.,0.); 
  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
    fParentPDGCode[npar] = -1; 
    fParentPythiaLine[npar] = -1;
  }
  for (Int_t i = 0; i < 4; i++){
    fQuarkPDGCode[i] = -1; 
    fQuarkPythiaLine[i] = -1; 
  }
  FillFromESD(muonTrack);
}

//============================================
AliMUONTrackLight::~AliMUONTrackLight()
{
/// Destructor
} 

//============================================
AliMUONTrackLight& AliMUONTrackLight::operator=(const AliMUONTrackLight& muonCopy)
{
  // check assignment to self
  if (this == &muonCopy) return *this;

  // base class assignment
  TObject::operator=(muonCopy);

  // assignment operator
  fPrec = muonCopy.fPrec; 
  fIsTriggered = muonCopy.fIsTriggered;
  fCharge = muonCopy.fCharge; 
  fChi2 = muonCopy.fChi2; 
  fCentr = muonCopy.fCentr;
  fPgen = muonCopy.fPgen; 
  fTrackPythiaLine = muonCopy.fTrackPythiaLine;
  fTrackPDGCode = muonCopy.fTrackPDGCode;
  fOscillation = muonCopy.fOscillation; 
  fNParents = muonCopy.fNParents;
  fWeight = muonCopy.fWeight;
  
  for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i]; 
  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
    fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar]; 
    fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
  }
  for (Int_t i = 0; i < 4; i++){
    fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i]; 
    fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i]; 
  }

  return *this;
}    

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

void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
  /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
  AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
  if (!trPar) {
    AliError("The track must contain the parameters at vertex");
    return;
  }
  this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
  this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz()); 
  this->SetTriggered(trackReco->GetMatchTrigger()); 
  
  Double_t xyz[3] = { trPar->GetNonBendingCoor(), 
		      trPar->GetBendingCoor(),
		      trPar->GetZ()};
  if (zvert!=-9999) xyz[2] = zvert;
  this->SetVertex(xyz); 
}

//============================================
void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
  /// computes prec and charge from ESD track
  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass(); 
  Double_t thetaX = muonTrack->GetThetaX();
  Double_t thetaY = muonTrack->GetThetaY();
  Double_t tanthx = TMath::Tan(thetaX);
  Double_t tanthy = TMath::Tan(thetaY);
  Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
  Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
  Double_t px     = pz * tanthx;
  Double_t py     = pz * tanthy;
  fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
  Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
  fPrec.SetPxPyPzE(px,py,pz,energy);
  // get the position
  fXYZ[0] = muonTrack->GetNonBendingCoor();
  fXYZ[1] = muonTrack->GetBendingCoor();
  if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
  else fXYZ[2] = zvert;
  // get the chi2 per d.o.f.
  fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
  fIsTriggered = muonTrack->GetMatchTrigger();
}

//============================================
void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){ 
  /// set the reconstructed 4-momentum, assuming the particle is a muon
  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass(); 
  Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
  fPrec.SetPxPyPzE(px,py,pz,energy);
}

//============================================
void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
  /// scans the muon history to determine parents pdg code and pythia line
  Int_t countP = -1;
  Int_t parents[10], parLine[10];
  Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother

  TParticle *mother;
  Int_t status=-1, pdg=-1;
  while(lineM >= 0){
    
    mother = stack->Particle(lineM); //direct mother of rec. track
    pdg = mother->GetPdgCode();//store PDG code of first mother
    // break if a string, gluon, quark or diquark is found 
    if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
    parents[++countP] = pdg;
    parLine[countP] = lineM;
    status = mother->GetStatusCode();//get its status code to check if oscillation occured
    if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
    lineM = mother->GetFirstMother();
  }
  //store all the fragmented parents in an array:
  for(int i = 0; i <= countP; i++){
    this->SetParentPDGCode(i,parents[countP-i]);
    this->SetParentPythiaLine(i,parLine[countP-i]);
  }
  fNParents = countP+1;
  countP = -1;

  //and store the lines of the string and further quarks in another array:
  while(lineM >= 0){
    mother = stack->Particle(lineM);
    pdg = mother->GetPdgCode();
    //now, get information before the fragmentation
    this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
    this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
    lineM = mother->GetFirstMother();
  }
  
  //check if in case of HF production, the string points to the correct end
  //and correct it in case of need:
  countP = 1;
  for(int par = 0; par < 4; par++){
    if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
      countP = par; //get the quark just before hadronisation
      break;
    }
  }
  if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
    if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){

      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
          this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
        );
      
      pdg = this->GetQuarkPDGCode(countP);
      Int_t line = this->GetQuarkPythiaLine(countP);
      this->ResetQuarkInfo();
      while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
                                                        //must coincide with the flavour of the last fragmented mother

	pdg = stack->Particle(++line)->GetPdgCode();    
      }
      //now, we have the correct string end and the correct line
      //continue to fill again all parents and corresponding lines
      while(line >= 0){
	mother = stack->Particle(line);//get again the mother
	pdg = mother->GetPdgCode();
	this->SetQuarkPythiaLine(countP, line);
	this->SetQuarkPDGCode(countP++, pdg);
	line = mother->GetFirstMother();
      }
      this->PrintInfo("h");
    }//mismatch
  }
}

//====================================
void AliMUONTrackLight::ResetQuarkInfo(){
  /// resets parton information
  for(int pos = 1; pos < 4; pos++){//[0] is the string
    this->SetQuarkPDGCode(pos,-1);
    this->SetQuarkPythiaLine(pos,-1);
  }
}

//====================================
Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
  /// checks if the particle is a B0 
  Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
  Bool_t answer = kFALSE;
  for(int i = 0; i < 2; i++){
    if(TMath::Abs(intTest) == bMes0[i]){
      answer = kTRUE;
      break;
    }
  }
  return answer;
}
//====================================
Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
  /// checks if mother is a resonance
  Int_t intTest = GetParentPDGCode(index); 
  // the resonance pdg code is built this way
  // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour 
  Int_t id=intTest%100000; 
  return (!((id-id%10)%110));
}
//====================================
Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
  /// returns the flavour of parent idParent (idParent=0 is the oldest 
  /// hadronized parent)
  Int_t pdg = GetParentPDGCode(idParent); 
  Int_t quark = TMath::Abs(pdg/100);
  if(quark > 9) quark = quark/10;
  return quark;
}

//====================================
void AliMUONTrackLight::PrintInfo(const Option_t* opt){
  /// prints information about the track: 
  /// - "H" muon's decay history
  /// - "K" muon kinematics
  /// - "A" all variables
  TString options(opt);
  options.ToUpper();

  if(options.Contains("H") || options.Contains("A")){ //muon decay history
    char *name= new char[100];
    TString pdg = "", line = "";
    for(int i = 3; i >= 0; i--){
      if(this->GetQuarkPythiaLine(i)>= 0){
	snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
	line += name;
	snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
	pdg += name;
      }
    }
    for(int i = 0; i < fNParents; i++){ 
      if(this->GetParentPythiaLine(i)>= 0){
	snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
	line += name;
	snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
	pdg += name;
      }
    }
    snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
    snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;

    printf("\nmuon's decay history:\n");
    printf(" PDG: %s\n", pdg.Data());
    printf("line: %s\n", line.Data());
  }
  if(options.Contains("K") || options.Contains("A")){ //muon kinematic

    Int_t charge = this->GetCharge();
    Double_t *vtx = this->GetVertex();
    TLorentzVector momRec = this->GetPRec();
    TLorentzVector momGen = this->GetPGen();
    printf("the track's charge is %d\n", charge);
    printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
    printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
    printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
    printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n", 
	   momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(), 
	   momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
  }
}
//====================================
Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
  /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
  Int_t pdg = this->GetParentPDGCode(idparent); 
  if (TMath::Abs(pdg)==211 || //pi+
      TMath::Abs(pdg)==321 || //K+
      TMath::Abs(pdg)==213 || //rho+
      TMath::Abs(pdg)==311 || //K0
      TMath::Abs(pdg)==313 || //K*0
      TMath::Abs(pdg)==323    //K*+
      ) { 
    return kTRUE;
  }
  else return kFALSE;
}
//====================================
Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
  /// check if the provided pdg code corresponds to a diquark 
  pdg = TMath::Abs(pdg);
  if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
  else return kFALSE;
}
Commit	Line	Data
55fd51b0	1	/**************************************************************************
55fd51b0	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
b88403f3	3	* SigmaEffect_thetadegrees *
55fd51b0	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 purpeateose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
3d1463c8	18	//-----------------------------------------------------------------------------
55fd51b0	19	// Compact information for the muon generated tracks in the MUON arm
	20	// useful at the last stage of the analysis chain
	21	// provides a link between the reconstructed track and the generated particle
	22	// stores kinematical information at gen. and rec. level and
	23	// the decay history of the muon, allowing the identification of the
	24	// mother process
	25	//
	26	// To be used together with AliMUONPairLight
3d1463c8	27	//
	28	// This class was prepared by INFN Cagliari, July 2006
	29	// (authors: H.Woehri, A.de Falco)
	30	//-----------------------------------------------------------------------------
55fd51b0	31
9bdda5f6	32	// 13 Nov 2007:
	33	// Added a temporary fix to FindRefTrack to be able to handle reconstructed tracks
	34	// generated from ESD muon track information. The problem is that the ESD data at
	35	// the moment only contains the first hit on chamber 1. Hopefully in the near future
	36	// this will be fixed and all hit information will be available.
	37	// - Artur Szostak <artursz@iafrica.com>
	38
55fd51b0	39	#include "AliMUONTrackLight.h"
	40	#include "AliMUONTrack.h"
	41	#include "AliMUONConstants.h"
bd5eb041	42	#include "AliMUONVTrackStore.h"
933daf4c	43	#include "AliMUONTrackParam.h"
55fd51b0	44
55fd51b0	45	#include "AliESDMuonTrack.h"
55fd51b0	46	#include "AliStack.h"
f202486b	47	#include "AliLog.h"
55fd51b0	48
55fd51b0	49	#include "TDatabasePDG.h"
55fd51b0	50	#include "TParticle.h"
	51	#include "TString.h"
	52
b785c4bd	53	#include <cstdio>
b785c4bd	54
55fd51b0	55	ClassImp(AliMUONTrackLight)
	56
	57	//===================================================================
	58
	59	AliMUONTrackLight::AliMUONTrackLight()
	60	: TObject(),
	61	fPrec(),
	62	fIsTriggered(kFALSE),
	63	fCharge(-999),
b88403f3	64	fChi2(-1),
55fd51b0	65	fCentr(-1),
	66	fPgen(),
	67	fTrackPythiaLine(-999),
	68	fTrackPDGCode(-999),
	69	fOscillation(kFALSE),
b88403f3	70	fNParents(0),
b88403f3	71	fWeight(1)
55fd51b0	72	{
	73	/// default constructor
	74	fPgen.SetPxPyPzE(0.,0.,0.,0.);
	75	fPrec.SetPxPyPzE(0.,0.,0.,0.);
	76	for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
	77	for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
	78	fParentPDGCode[npar] = -1;
	79	fParentPythiaLine[npar] = -1;
	80	}
	81	for (Int_t i = 0; i < 4; i++){
	82	fQuarkPDGCode[i] = -1;
	83	fQuarkPythiaLine[i] = -1;
	84	}
	85	}
	86
	87	//============================================
	88	AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
	89	: TObject(muonCopy),
	90	fPrec(muonCopy.fPrec),
	91	fIsTriggered(muonCopy.fIsTriggered),
	92	fCharge(muonCopy.fCharge),
b88403f3	93	fChi2(muonCopy.fChi2),
55fd51b0	94	fCentr(muonCopy.fCentr),
	95	fPgen(muonCopy.fPgen),
	96	fTrackPythiaLine(muonCopy.fTrackPythiaLine),
	97	fTrackPDGCode(muonCopy.fTrackPDGCode),
	98	fOscillation(muonCopy.fOscillation),
b88403f3	99	fNParents(muonCopy.fNParents),
b88403f3	100	fWeight(muonCopy.fWeight)
55fd51b0	101	{
	102	/// copy constructor
	103	for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
	104	for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
	105	fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
	106	fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
	107	}
	108	for (Int_t i = 0; i < 4; i++){
	109	fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
	110	fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
	111	}
	112	}
	113
	114	//============================================
	115	AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
	116	: TObject(),
	117	fPrec(),
	118	fIsTriggered(kFALSE),
	119	fCharge(-999),
b88403f3	120	fChi2(-1),
55fd51b0	121	fCentr(-1),
	122	fPgen(),
	123	fTrackPythiaLine(-999),
	124	fTrackPDGCode(-999),
	125	fOscillation(kFALSE),
b88403f3	126	fNParents(0),
b88403f3	127	fWeight(1)
55fd51b0	128	{
55fd51b0	129	/// constructor
f202486b	130	fPgen.SetPxPyPzE(0.,0.,0.,0.);
	131	for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
	132	fParentPDGCode[npar] = -1;
	133	fParentPythiaLine[npar] = -1;
	134	}
	135	for (Int_t i = 0; i < 4; i++){
	136	fQuarkPDGCode[i] = -1;
	137	fQuarkPythiaLine[i] = -1;
	138	}
	139	FillFromESD(muonTrack);
55fd51b0	140	}
	141
	142	//============================================
b88403f3	143	AliMUONTrackLight::~AliMUONTrackLight()
	144	{
	145	/// Destructor
	146	}
	147
e733e3dd	148	//============================================
	149	AliMUONTrackLight& AliMUONTrackLight::operator=(const AliMUONTrackLight& muonCopy)
	150	{
	151	// check assignment to self
	152	if (this == &muonCopy) return *this;
	153
	154	// base class assignment
	155	TObject::operator=(muonCopy);
	156
	157	// assignment operator
	158	fPrec = muonCopy.fPrec;
	159	fIsTriggered = muonCopy.fIsTriggered;
	160	fCharge = muonCopy.fCharge;
	161	fChi2 = muonCopy.fChi2;
	162	fCentr = muonCopy.fCentr;
	163	fPgen = muonCopy.fPgen;
	164	fTrackPythiaLine = muonCopy.fTrackPythiaLine;
	165	fTrackPDGCode = muonCopy.fTrackPDGCode;
	166	fOscillation = muonCopy.fOscillation;
	167	fNParents = muonCopy.fNParents;
	168	fWeight = muonCopy.fWeight;
	169
	170	for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
	171	for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
	172	fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
	173	fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
	174	}
	175	for (Int_t i = 0; i < 4; i++){
	176	fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
	177	fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
	178	}
72421d27	179
72421d27	180	return *this;
e733e3dd	181	}
e733e3dd	182
b88403f3	183	//============================================
	184
	185	void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
	186	/// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
f202486b	187	AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
	188	if (!trPar) {
	189	AliError("The track must contain the parameters at vertex");
	190	return;
61fed964	191	}
f202486b	192	this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
f202486b	193	this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
b88403f3	194	this->SetTriggered(trackReco->GetMatchTrigger());
b88403f3	195
f202486b	196	Double_t xyz[3] = { trPar->GetNonBendingCoor(),
	197	trPar->GetBendingCoor(),
	198	trPar->GetZ()};
61fed964	199	if (zvert!=-9999) xyz[2] = zvert;
b88403f3	200	this->SetVertex(xyz);
	201	}
	202
	203	//============================================
	204	void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
55fd51b0	205	/// computes prec and charge from ESD track
	206	Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
	207	Double_t thetaX = muonTrack->GetThetaX();
	208	Double_t thetaY = muonTrack->GetThetaY();
	209	Double_t tanthx = TMath::Tan(thetaX);
	210	Double_t tanthy = TMath::Tan(thetaY);
	211	Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
	212	Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
	213	Double_t px = pz * tanthx;
	214	Double_t py = pz * tanthy;
	215	fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
	216	Double_t energy = TMath::Sqrt(mumass * mumass + pxpx + pypy + pz*pz);
	217	fPrec.SetPxPyPzE(px,py,pz,energy);
b88403f3	218	// get the position
	219	fXYZ[0] = muonTrack->GetNonBendingCoor();
	220	fXYZ[1] = muonTrack->GetBendingCoor();
	221	if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
	222	else fXYZ[2] = zvert;
	223	// get the chi2 per d.o.f.
	224	fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
	225	fIsTriggered = muonTrack->GetMatchTrigger();
55fd51b0	226	}
	227
	228	//============================================
	229	void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
	230	/// set the reconstructed 4-momentum, assuming the particle is a muon
	231	Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
	232	Double_t energy = TMath::Sqrt(mumass * mumass + pxpx + pypy + pz*pz);
	233	fPrec.SetPxPyPzE(px,py,pz,energy);
	234	}
	235
	236	//============================================
b88403f3	237	void AliMUONTrackLight::FillMuonHistory(AliStack stack, TParticle part){
55fd51b0	238	/// scans the muon history to determine parents pdg code and pythia line
	239	Int_t countP = -1;
	240	Int_t parents[10], parLine[10];
	241	Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
	242
55fd51b0	243	TParticle *mother;
	244	Int_t status=-1, pdg=-1;
	245	while(lineM >= 0){
b88403f3	246
55fd51b0	247	mother = stack->Particle(lineM); //direct mother of rec. track
55fd51b0	248	pdg = mother->GetPdgCode();//store PDG code of first mother
b88403f3	249	// break if a string, gluon, quark or diquark is found
b88403f3	250	if(pdg == 92 \|\| pdg == 21 \|\| TMath::Abs(pdg) < 10 \|\| IsDiquark(pdg)) break;
55fd51b0	251	parents[++countP] = pdg;
	252	parLine[countP] = lineM;
	253	status = mother->GetStatusCode();//get its status code to check if oscillation occured
	254	if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
	255	lineM = mother->GetFirstMother();
	256	}
	257	//store all the fragmented parents in an array:
	258	for(int i = 0; i <= countP; i++){
	259	this->SetParentPDGCode(i,parents[countP-i]);
	260	this->SetParentPythiaLine(i,parLine[countP-i]);
	261	}
55fd51b0	262	fNParents = countP+1;
55fd51b0	263	countP = -1;
b88403f3	264
55fd51b0	265	//and store the lines of the string and further quarks in another array:
	266	while(lineM >= 0){
	267	mother = stack->Particle(lineM);
	268	pdg = mother->GetPdgCode();
	269	//now, get information before the fragmentation
	270	this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
	271	this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
	272	lineM = mother->GetFirstMother();
	273	}
b88403f3	274
55fd51b0	275	//check if in case of HF production, the string points to the correct end
	276	//and correct it in case of need:
	277	countP = 1;
b88403f3	278	for(int par = 0; par < 4; par++){
	279	if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
	280	countP = par; //get the quark just before hadronisation
	281	break;
	282	}
	283	}
55fd51b0	284	if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 \|\| this->GetParentFlavour(0)==5)){
	285	if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
	286
9bdda5f6	287	AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
	288	this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
	289	);
	290
7cdc832b	291	pdg = this->GetQuarkPDGCode(countP);
7cdc832b	292	Int_t line = this->GetQuarkPythiaLine(countP);
55fd51b0	293	this->ResetQuarkInfo();
	294	while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
	295	//must coincide with the flavour of the last fragmented mother
	296
	297	pdg = stack->Particle(++line)->GetPdgCode();
	298	}
	299	//now, we have the correct string end and the correct line
	300	//continue to fill again all parents and corresponding lines
	301	while(line >= 0){
	302	mother = stack->Particle(line);//get again the mother
	303	pdg = mother->GetPdgCode();
	304	this->SetQuarkPythiaLine(countP, line);
	305	this->SetQuarkPDGCode(countP++, pdg);
	306	line = mother->GetFirstMother();
	307	}
b88403f3	308	this->PrintInfo("h");
b88403f3	309	}//mismatch
55fd51b0	310	}
55fd51b0	311	}
b88403f3	312
55fd51b0	313	//====================================
	314	void AliMUONTrackLight::ResetQuarkInfo(){
	315	/// resets parton information
	316	for(int pos = 1; pos < 4; pos++){//[0] is the string
	317	this->SetQuarkPDGCode(pos,-1);
	318	this->SetQuarkPythiaLine(pos,-1);
	319	}
	320	}
	321
	322	//====================================
	323	Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
	324	/// checks if the particle is a B0
	325	Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
	326	Bool_t answer = kFALSE;
	327	for(int i = 0; i < 2; i++){
	328	if(TMath::Abs(intTest) == bMes0[i]){
	329	answer = kTRUE;
	330	break;
	331	}
	332	}
	333	return answer;
	334	}
	335	//====================================
	336	Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
	337	/// checks if mother is a resonance
	338	Int_t intTest = GetParentPDGCode(index);
	339	// the resonance pdg code is built this way
	340	// x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
	341	Int_t id=intTest%100000;
	342	return (!((id-id%10)%110));
	343	}
	344	//====================================
	345	Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
	346	/// returns the flavour of parent idParent (idParent=0 is the oldest
	347	/// hadronized parent)
	348	Int_t pdg = GetParentPDGCode(idParent);
	349	Int_t quark = TMath::Abs(pdg/100);
	350	if(quark > 9) quark = quark/10;
	351	return quark;
	352	}
	353
	354	//====================================
57e2ad1a	355	void AliMUONTrackLight::PrintInfo(const Option_t* opt){
55fd51b0	356	/// prints information about the track:
	357	/// - "H" muon's decay history
	358	/// - "K" muon kinematics
	359	/// - "A" all variables
	360	TString options(opt);
	361	options.ToUpper();
	362
	363	if(options.Contains("H") \|\| options.Contains("A")){ //muon decay history
	364	char *name= new char[100];
	365	TString pdg = "", line = "";
	366	for(int i = 3; i >= 0; i--){
	367	if(this->GetQuarkPythiaLine(i)>= 0){
b785c4bd	368	snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
55fd51b0	369	line += name;
b785c4bd	370	snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
55fd51b0	371	pdg += name;
	372	}
	373	}
	374	for(int i = 0; i < fNParents; i++){
	375	if(this->GetParentPythiaLine(i)>= 0){
b785c4bd	376	snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
55fd51b0	377	line += name;
b785c4bd	378	snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
55fd51b0	379	pdg += name;
	380	}
	381	}
b785c4bd	382	snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
b785c4bd	383	snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
55fd51b0	384
	385	printf("\nmuon's decay history:\n");
	386	printf(" PDG: %s\n", pdg.Data());
	387	printf("line: %s\n", line.Data());
	388	}
	389	if(options.Contains("K") \|\| options.Contains("A")){ //muon kinematic
	390
	391	Int_t charge = this->GetCharge();
	392	Double_t *vtx = this->GetVertex();
	393	TLorentzVector momRec = this->GetPRec();
	394	TLorentzVector momGen = this->GetPGen();
	395	printf("the track's charge is %d\n", charge);
	396	printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
	397	printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
	398	printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
	399	printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
	400	momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
	401	momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
	402	}
	403	}
b88403f3	404	//====================================
	405	Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
	406	/// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
	407	Int_t pdg = this->GetParentPDGCode(idparent);
	408	if (TMath::Abs(pdg)==211 \|\| //pi+
	409	TMath::Abs(pdg)==321 \|\| //K+
	410	TMath::Abs(pdg)==213 \|\| //rho+
	411	TMath::Abs(pdg)==311 \|\| //K0
	412	TMath::Abs(pdg)==313 \|\| //K*0
	413	TMath::Abs(pdg)==323 //K*+
	414	) {
	415	return kTRUE;
	416	}
	417	else return kFALSE;
	418	}
	419	//====================================
9ee1d6ff	420	Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
b88403f3	421	/// check if the provided pdg code corresponds to a diquark
	422	pdg = TMath::Abs(pdg);
	423	if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
	424	else return kFALSE;
	425	}