[u/mrichter/AliRoot.git] / MUON / AliMUONPairLight.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$ */

//-----------------------------------------------------------------------------
//This class was prepared by INFN Cagliari, July 2006
//(authors: H.Woehri, A.de Falco)
// 
// Compact information for the generated muon pairs in the MUON arm 
// useful at the last stage of the analysis chain
// Pairs are built with two AliMUONTrackLight objects 
// Using the class AliMUONTrackLight this class combines the decay
// information ("history") of the reconstructed tracks and fills
// a series of flags for the formed reconstructed dimuon:
// fIsCorrelated, fCreationProcess, fIsFeedDown, ...
// for information about the dimuon, use PrintInfo with the appropriate
// printflag
// To be used together with AliMUONTrackLight
//-----------------------------------------------------------------------------


//MUON classes
#include "AliMUONPairLight.h"
//Root classes
#include "TString.h"

ClassImp(AliMUONPairLight) 

//====================================
AliMUONPairLight::AliMUONPairLight() : 
  TObject(), 
  fMu0(),
  fMu1(), 
  fCreationProcess(-999),
  fIsCorrelated(kFALSE), 
  fCauseOfCorrelation (-1),
  fIsFeedDown(kFALSE)
{
  /// default constructor
  ; 
}

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

AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy) 
  : TObject(dimuCopy),
    fMu0(dimuCopy.fMu0),
    fMu1(dimuCopy.fMu1), 
    fCreationProcess(dimuCopy.fCreationProcess),
    fIsCorrelated(dimuCopy.fIsCorrelated), 
    fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
    fIsFeedDown(dimuCopy.fIsFeedDown)
{ 
/// copy constructor
///   fMu0 = AliMUONTrackLight(dimuCopy.fMu0); 
///   fMu1 = AliMUONTrackLight(dimuCopy.fMu1); 
///   fIsCorrelated = dimuCopy.fIsCorrelated;
///   fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
///   fCreationProcess = dimuCopy.fCreationProcess;
///   fIsFeedDown = dimuCopy.fIsFeedDown;
  ;
}

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

AliMUONPairLight::~AliMUONPairLight(){
  /// destructor
}

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

Bool_t AliMUONPairLight::IsAResonance(){
  /// checks if muon pair comes from a resonance decay  
  if (!fIsCorrelated) return kFALSE;   //if muons not correlated, cannot be a resonance
  //if muons are correlated, check if the PDG of the
  //common mother is a resonance
  Int_t nparents0 = fMu0.GetNParents(); 
  Int_t nparents1 = fMu1.GetNParents(); 

  Int_t minP = TMath::Min(nparents0, nparents1);
  for (Int_t i = 0 ; i < minP; i++) { 
    if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) && 
	fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
      if (nparents0-1-i) SetFeedDown(nparents0-1-i);
      return kTRUE;
    }
  }
  return kFALSE; 
}

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

AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index)  { 
  /// return muon 0 or 1
   if (index==0) return &fMu0;
   else if (index==1) return &fMu1; 
   else{ printf ("Index can be either 0 or 1\n"); return 0;}
   //   else return &fMu1; 
}

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

Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) { 
  /// return muon mother pdg code
  if (imuon==0) return fMu0.GetParentPDGCode(mother); 
  else if (imuon==1) return fMu1.GetParentPDGCode(mother); 
  else { printf ("Index must be only 0 or 1\n"); return -999; } 
}

//====================================
void AliMUONPairLight::SetProcess(){
  /// finds the process related to the muon pair (open charm/beauty, resonance, 
  /// uncorrelated...) 

  AliMUONTrackLight *mu1 = &fMu0;
  AliMUONTrackLight *mu2 = &fMu1;

  // check if the two muons are correlated
  // first check if they come from the same hadron (resonance or beauty/charm meson)
  Int_t npar1 = mu1->GetNParents(); 
  Int_t npar2 = mu2->GetNParents(); 
  for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) { 
    Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
    for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) { 
      Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
      if(lineMo1 == lineMo2) { 
	//reject "diquark" mothers
	if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
// 	if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
	this->SetCorrelated(kTRUE); 
	this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
	if(!IsAResonance()) fCreationProcess = 3; 
	else fCreationProcess = -1;
	return;
      }
    }
  }

  //now, check if we have a correlated pi/K:
  if(this->IsDimuonFromCorrPiK()){
    this->SetCorrelated(kTRUE); 
    this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
    fCreationProcess = -1;
  }

  // if Open Beauty/Charm we can have 3 creation processes 
  // (pair creation [0], gluon splitting [1] or flavour excitation [2])
  // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same 
  Int_t flavPar1 = mu1->GetParentFlavour(0);
  Int_t flavPar2 = mu2->GetParentFlavour(0);
  for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) { 
    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
    for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) { 
      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
      if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
	//now, check also that the string fragmented into two hadrons
	//of the same flavour (string usually splits into many hadrons
	//among which there are mostly soft particles)
	if(flavPar1 == flavPar2){
	  this->SetCorrelated(kTRUE); 
	  if(GetCauseOfCorrelation() == -1)
	    this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));

	  fCreationProcess = 1; 
	  return;
	}
      }
    }
  }

  Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
  Int_t line2 = mu2->GetQuarkPythiaLine(2); 

  Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
  Int_t flavourLine6or7[2] = {-1, -1};
  // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
  // are filled with a Q and Qbar
  for (Int_t imoth1 = 3; imoth1>=0; imoth1--) { 
    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
    Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
    if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
      line6or7[0] = imoth1;
      flavourLine6or7[0] = flavour1;
    }
    for (Int_t imoth2 = 3; imoth2>=0; imoth2--) { 
      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
      Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
      if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
	line6or7[1] = imoth2;
	flavourLine6or7[1] = flavour2;
      }
      if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
	 (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
	 (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
	 (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour 
	this->SetCorrelated(kTRUE);
	fCreationProcess = 0; 
	return;
      }
    }
  }

  // 3.)flavour excitation: if pythia line 6 of one track *and* pythia line 7 of second track
  // are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
  Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
  Int_t flavourLine2or3[2] = {-1, -1};
  for (Int_t imoth1 = 3; imoth1>=0; imoth1--) { 
    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
    Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
    if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
      line2or3[0] = imoth1;
      flavourLine2or3[0] = flavour1;
    }
    for (Int_t imoth2 = 3; imoth2>=0; imoth2--) { 
      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
      Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
      if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
	line2or3[1] = imoth2;
	flavourLine2or3[1] = flavour2;
      }
      if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4  || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
	  (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
	 ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) &&  //or the same,
	  (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
	//now, check also that the string fragmented into two hadrons
	//of the same flavour (string usually splits into many hadrons
	//among which there are mostly soft particles)
	if(flavPar1 == flavPar2){
	  this->SetCorrelated(kTRUE);
	  fCreationProcess = 2;
	  return;
	}
      }
    }
  } 

  //now flag (rare) processes in which only the incoming parton in line 2 or 3
  //radiates a gluon which produces a QQbar pair:
  //exclude the light quarks
  if(line1 == line2 && (line1 == 2 || line1 == 3)){
    if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
       (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){

      //now, check also that the string fragmented into two hadrons
      //of the same flavour (string usually splits into many hadrons
      //among which there are mostly soft particles)
      if(flavPar1 == flavPar2){

	this->SetCorrelated(kTRUE);
	fCreationProcess = 1;
	if(GetCauseOfCorrelation() == -1){
	  this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
	}
	return;
      }
    }
  }

  //in initial-state-radiation produced QQbar events the "mother quark"
  //is acknowledged as the second quark [1] and sits in line 2 or 3
  //is part of gluon splitting
  line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
  line2 = mu2->GetQuarkPythiaLine(1); 
  if(line1 == line2 && (line1 == 2 || line1 == 3)){
    if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
       (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){

      //now, check also that the string fragmented into two hadrons
      //of the same flavour (string usually splits into many hadrons
      //among which there are mostly soft particles)
      if(flavPar1 == flavPar2){
	
	this->SetCorrelated(kTRUE);
	fCreationProcess = 1;
	if(GetCauseOfCorrelation() == -1){
	  this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
	}
	return;
      }
    }
  }

  //in final-state-radiation produced QQbar events the "mother quark"
  //is acknowledged as the first quark [1] and sits in line 6 or 7
  //is part of gluon splitting
  line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark 
  line2 = mu2->GetQuarkPythiaLine(1); 
  if(line1 == line2 && (line1 == 6 || line1 == 7)){
    if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
       (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){

      //now, check also that the string fragmented into two hadrons
      //of the same flavour (string usually splits into many hadrons
      //among which there are mostly soft particles)
      if(flavPar1 == flavPar2){
	
	this->SetCorrelated(kTRUE);
	fCreationProcess = 1;
	if(GetCauseOfCorrelation() == -1){
	  this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
	}
	return;
      }
    }
  }
}

//====================================
void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
  /// set the two muons 
  fMu0 = mu0; 
  fMu1 = mu1; 
  this->SetProcess();
} 

//====================================
void AliMUONPairLight::PrintInfo(Option_t* opt){
  /// print information about muon pairs
  /// Options: 
  /// - "H" single muons' decay histories
  /// - "K" dimuon kinematics
  /// - "F" dimuon flags
  /// - "A" all variables
  TString options(opt);
  options.ToUpper();

  if(options.Contains("H") || options.Contains("A")){//muon decay histories

    AliMUONTrackLight *mu1 = &fMu0;
    AliMUONTrackLight *mu2 = &fMu1;

    printf("========= History =======================\n");
    printf("first muon");
    mu1->PrintInfo("H");
    printf("second muon");
    mu2->PrintInfo("H");
    printf("=========================================\n");
  }
  if(options.Contains("F") || options.Contains("A")){//flags
    printf("the flags set for this muon pair are:\n");
    printf("=====================================\n");
    if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
    fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
    if(IsOpenCharm()) printf("(*) correlated open charm: ");
    if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
    if(IsOpenCharm() || IsOpenBeauty()){
      switch(fCreationProcess){
      case 0:
	printf("pair creation");
	break;
      case 1:
	printf("gluon splitting");
	break;
      case 2:
	printf("flavour excitation");
	break;
      case 3:
	printf("both muons come from same fragmented mother");
	break;
      }
      if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation()) 
	printf("... where oscillation occured\n");
      else{
	if(IsOpenBeauty())
	  printf(" (no oscillation)\n");
	else
	  printf("\n");
      }
    }
    IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
    fIsFeedDown ? printf("(*) mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(0,fMu0.GetNParents()-2), this->GetMuonMotherPDG(0,fMu0.GetNParents()-1)) : printf("(*) no feed-down\n");
    printf("=====================================\n");
  }
  if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
    Double_t *vtx = this->GetMuon(0)->GetVertex();
    TLorentzVector momRec = this->GetPRec();
    TLorentzVector momGen = this->GetPGen();
    printf("the dimuon charge is %d\n", this->GetCharge());
    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());
    //rapidity, pT, angles, ...
    printf("Rec. variables: mass %1.3f, pT %1.3f, pseudo-rapidity %1.3f, openingAngle %1.3f (%1.3f degree), theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n", 
	   momRec.M(), momRec.Pt(), momRec.Eta(), 
	   TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(), 
	   momRec.Theta(), 180./TMath::Pi() * momRec.Theta(), 
	   momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
  }
}

//====================================
Double_t AliMUONPairLight::GetOpeningAngle() { 
  /// opening angle between the two muons in the lab frame (in degrees)
  TLorentzVector pRecMu0 =  fMu0.GetPRec();
  TLorentzVector pRecMu1 =  fMu1.GetPRec();
  TVector3 pRecMu03 = pRecMu0.Vect();
  TVector3 pRecMu13 = pRecMu1.Vect();
  Double_t scalar = pRecMu03.Dot(pRecMu13);
  Double_t modMu0 = pRecMu03.Mag();
  Double_t modMu1 = pRecMu13.Mag();
  Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
  return theta; 
}
//================================================
Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
  ///check if we have a correlated pi/K

  AliMUONTrackLight *mu0 = this->GetMuon(0), *mu1 = this->GetMuon(1);
  Bool_t fromSameLine = kFALSE;
  if (mu0->IsParentPionOrKaon() &&
      mu1->IsParentPionOrKaon() &&
      mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
      ) fromSameLine = kTRUE;

  return fromSameLine;
}
Commit	Line	Data
55fd51b0	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* SigmaEffect_thetadegrees *
	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	//This class was prepared by INFN Cagliari, July 2006
	20	//(authors: H.Woehri, A.de Falco)
	21	//
	22	// Compact information for the generated muon pairs in the MUON arm
	23	// useful at the last stage of the analysis chain
	24	// Pairs are built with two AliMUONTrackLight objects
	25	// Using the class AliMUONTrackLight this class combines the decay
	26	// information ("history") of the reconstructed tracks and fills
	27	// a series of flags for the formed reconstructed dimuon:
	28	// fIsCorrelated, fCreationProcess, fIsFeedDown, ...
	29	// for information about the dimuon, use PrintInfo with the appropriate
	30	// printflag
	31	// To be used together with AliMUONTrackLight
3d1463c8	32	//-----------------------------------------------------------------------------
55fd51b0	33
	34
	35	//MUON classes
	36	#include "AliMUONPairLight.h"
	37	//Root classes
	38	#include "TString.h"
	39
	40	ClassImp(AliMUONPairLight)
	41
	42	//====================================
	43	AliMUONPairLight::AliMUONPairLight() :
	44	TObject(),
	45	fMu0(),
	46	fMu1(),
	47	fCreationProcess(-999),
	48	fIsCorrelated(kFALSE),
	49	fCauseOfCorrelation (-1),
	50	fIsFeedDown(kFALSE)
	51	{
	52	/// default constructor
	53	;
	54	}
	55
	56	//====================================
	57
	58	AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy)
	59	: TObject(dimuCopy),
	60	fMu0(dimuCopy.fMu0),
0060b2f7	61	fMu1(dimuCopy.fMu1),
55fd51b0	62	fCreationProcess(dimuCopy.fCreationProcess),
	63	fIsCorrelated(dimuCopy.fIsCorrelated),
	64	fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
	65	fIsFeedDown(dimuCopy.fIsFeedDown)
	66	{
	67	/// copy constructor
	68	/// fMu0 = AliMUONTrackLight(dimuCopy.fMu0);
	69	/// fMu1 = AliMUONTrackLight(dimuCopy.fMu1);
	70	/// fIsCorrelated = dimuCopy.fIsCorrelated;
	71	/// fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
	72	/// fCreationProcess = dimuCopy.fCreationProcess;
	73	/// fIsFeedDown = dimuCopy.fIsFeedDown;
	74	;
	75	}
	76
	77	//====================================
	78
	79	AliMUONPairLight::~AliMUONPairLight(){
	80	/// destructor
	81	}
	82
	83	//====================================
	84
	85	Bool_t AliMUONPairLight::IsAResonance(){
	86	/// checks if muon pair comes from a resonance decay
	87	if (!fIsCorrelated) return kFALSE; //if muons not correlated, cannot be a resonance
	88	//if muons are correlated, check if the PDG of the
	89	//common mother is a resonance
b88403f3	90	Int_t nparents0 = fMu0.GetNParents();
b88403f3	91	Int_t nparents1 = fMu1.GetNParents();
55fd51b0	92
b88403f3	93	Int_t minP = TMath::Min(nparents0, nparents1);
	94	for (Int_t i = 0 ; i < minP; i++) {
	95	if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) &&
	96	fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
	97	if (nparents0-1-i) SetFeedDown(nparents0-1-i);
	98	return kTRUE;
	99	}
	100	}
	101	return kFALSE;
55fd51b0	102	}
	103
	104	//====================================
	105
	106	AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index) {
	107	/// return muon 0 or 1
	108	if (index==0) return &fMu0;
	109	else if (index==1) return &fMu1;
	110	else{ printf ("Index can be either 0 or 1\n"); return 0;}
	111	// else return &fMu1;
	112	}
	113
	114	//====================================
	115
	116	Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
	117	/// return muon mother pdg code
	118	if (imuon==0) return fMu0.GetParentPDGCode(mother);
	119	else if (imuon==1) return fMu1.GetParentPDGCode(mother);
	120	else { printf ("Index must be only 0 or 1\n"); return -999; }
	121	}
	122
	123	//====================================
	124	void AliMUONPairLight::SetProcess(){
	125	/// finds the process related to the muon pair (open charm/beauty, resonance,
	126	/// uncorrelated...)
b88403f3	127
55fd51b0	128	AliMUONTrackLight *mu1 = &fMu0;
	129	AliMUONTrackLight *mu2 = &fMu1;
	130
55fd51b0	131	// check if the two muons are correlated
	132	// first check if they come from the same hadron (resonance or beauty/charm meson)
	133	Int_t npar1 = mu1->GetNParents();
	134	Int_t npar2 = mu2->GetNParents();
55fd51b0	135	for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
55fd51b0	136	Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
55fd51b0	137	for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
	138	Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
	139	if(lineMo1 == lineMo2) {
b88403f3	140	//reject "diquark" mothers
	141	if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
	142	// if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
55fd51b0	143	this->SetCorrelated(kTRUE);
	144	this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
	145	if(!IsAResonance()) fCreationProcess = 3;
	146	else fCreationProcess = -1;
b88403f3	147	return;
55fd51b0	148	}
55fd51b0	149	}
91e620e2	150	}
	151
	152	//now, check if we have a correlated pi/K:
	153	if(this->IsDimuonFromCorrPiK()){
	154	this->SetCorrelated(kTRUE);
	155	this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
	156	fCreationProcess = -1;
	157	}
	158
55fd51b0	159	// if Open Beauty/Charm we can have 3 creation processes
55fd51b0	160	// (pair creation [0], gluon splitting [1] or flavour excitation [2])
55fd51b0	161	// 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
b88403f3	162	Int_t flavPar1 = mu1->GetParentFlavour(0);
	163	Int_t flavPar2 = mu2->GetParentFlavour(0);
	164	for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) {
	165	Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
	166	for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) {
	167	Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
	168	if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
	169	//now, check also that the string fragmented into two hadrons
	170	//of the same flavour (string usually splits into many hadrons
	171	//among which there are mostly soft particles)
	172	if(flavPar1 == flavPar2){
	173	this->SetCorrelated(kTRUE);
	174	if(GetCauseOfCorrelation() == -1)
	175	this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
	176
	177	fCreationProcess = 1;
	178	return;
	179	}
	180	}
55fd51b0	181	}
55fd51b0	182	}
91e620e2	183
55fd51b0	184	Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
55fd51b0	185	Int_t line2 = mu2->GetQuarkPythiaLine(2);
b88403f3	186
	187	Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
	188	Int_t flavourLine6or7[2] = {-1, -1};
	189	// 2.) pair creation: if pythia line 6 of one track and pythia line 7 of second track
	190	// are filled with a Q and Qbar
	191	for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
	192	Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
	193	Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
	194	if(lineMo1 == 6 \|\| lineMo1 == 7){ //track 0 has a mother in line 6 or 7
	195	line6or7[0] = imoth1;
	196	flavourLine6or7[0] = flavour1;
	197	}
	198	for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
	199	Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
	200	Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
	201	if(lineMo2 == 6 \|\| lineMo2 == 7){ //track 1 has a mother in line 6 or 7
	202	line6or7[1] = imoth2;
	203	flavourLine6or7[1] = flavour2;
	204	}
	205	if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
	206	(flavourLine6or7[0] == 4 \|\| flavourLine6or7[0] == 5) && //this entry must be a c or b quark
	207	(flavourLine6or7[1] == 4 \|\| flavourLine6or7[1] == 5) && // == " ==
	208	(flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour
	209	this->SetCorrelated(kTRUE);
	210	fCreationProcess = 0;
	211	return;
	212	}
55fd51b0	213	}
	214	}
	215
91e620e2	216	// 3.)flavour excitation: if pythia line 6 of one track and pythia line 7 of second track
91e620e2	217	// are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
b88403f3	218	Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
	219	Int_t flavourLine2or3[2] = {-1, -1};
	220	for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
	221	Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
	222	Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
	223	if(lineMo1 == 2 \|\| lineMo1 == 3){ //track 0 has a mother in line 2 or 3
	224	line2or3[0] = imoth1;
	225	flavourLine2or3[0] = flavour1;
	226	}
	227	for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
	228	Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
	229	Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
	230	if(lineMo2 == 2 \|\| lineMo2 == 3){ //track 1 has a mother in line 2 or 3
	231	line2or3[1] = imoth2;
	232	flavourLine2or3[1] = flavour2;
	233	}
	234	if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4 \|\| flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
	235	(line2or3[1] > 0 && (flavourLine2or3[1] == 21 \|\| flavourLine2or3[1] < 10))) \|\| //second track has a g/q in line 2 or 3
	236	((line6or7[1] > 0 && (flavourLine6or7[1] == 4 \|\| flavourLine6or7[1] == 5)) && //or the same,
	237	(line2or3[0] > 0 && (flavourLine2or3[0] == 21 \|\| flavourLine2or3[0] < 10)))){ // swapping the track's indices
	238	//now, check also that the string fragmented into two hadrons
	239	//of the same flavour (string usually splits into many hadrons
	240	//among which there are mostly soft particles)
	241	if(flavPar1 == flavPar2){
	242	this->SetCorrelated(kTRUE);
	243	fCreationProcess = 2;
	244	return;
	245	}
	246	}
	247	}
	248	}
	249
	250	//now flag (rare) processes in which only the incoming parton in line 2 or 3
	251	//radiates a gluon which produces a QQbar pair:
	252	//exclude the light quarks
	253	if(line1 == line2 && (line1 == 2 \|\| line1 == 3)){
	254	if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) \|\|
	255	(TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
	256
	257	//now, check also that the string fragmented into two hadrons
	258	//of the same flavour (string usually splits into many hadrons
	259	//among which there are mostly soft particles)
	260	if(flavPar1 == flavPar2){
	261
	262	this->SetCorrelated(kTRUE);
	263	fCreationProcess = 1;
	264	if(GetCauseOfCorrelation() == -1){
	265	this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
	266	}
	267	return;
	268	}
	269	}
55fd51b0	270	}
91e620e2	271
	272	//in initial-state-radiation produced QQbar events the "mother quark"
	273	//is acknowledged as the second quark [1] and sits in line 2 or 3
	274	//is part of gluon splitting
	275	line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
	276	line2 = mu2->GetQuarkPythiaLine(1);
	277	if(line1 == line2 && (line1 == 2 \|\| line1 == 3)){
	278	if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) \|\|
	279	(TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
	280
	281	//now, check also that the string fragmented into two hadrons
	282	//of the same flavour (string usually splits into many hadrons
	283	//among which there are mostly soft particles)
	284	if(flavPar1 == flavPar2){
	285
	286	this->SetCorrelated(kTRUE);
	287	fCreationProcess = 1;
	288	if(GetCauseOfCorrelation() == -1){
	289	this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
	290	}
	291	return;
	292	}
	293	}
	294	}
	295
	296	//in final-state-radiation produced QQbar events the "mother quark"
	297	//is acknowledged as the first quark [1] and sits in line 6 or 7
	298	//is part of gluon splitting
	299	line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark
	300	line2 = mu2->GetQuarkPythiaLine(1);
	301	if(line1 == line2 && (line1 == 6 \|\| line1 == 7)){
	302	if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) \|\|
	303	(TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
	304
	305	//now, check also that the string fragmented into two hadrons
	306	//of the same flavour (string usually splits into many hadrons
	307	//among which there are mostly soft particles)
	308	if(flavPar1 == flavPar2){
	309
	310	this->SetCorrelated(kTRUE);
	311	fCreationProcess = 1;
	312	if(GetCauseOfCorrelation() == -1){
	313	this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
	314	}
	315	return;
	316	}
	317	}
	318	}
55fd51b0	319	}
	320
	321	//====================================
	322	void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
	323	/// set the two muons
	324	fMu0 = mu0;
	325	fMu1 = mu1;
	326	this->SetProcess();
	327	}
	328
	329	//====================================
	330	void AliMUONPairLight::PrintInfo(Option_t* opt){
	331	/// print information about muon pairs
	332	/// Options:
	333	/// - "H" single muons' decay histories
	334	/// - "K" dimuon kinematics
	335	/// - "F" dimuon flags
	336	/// - "A" all variables
	337	TString options(opt);
	338	options.ToUpper();
	339
	340	if(options.Contains("H") \|\| options.Contains("A")){//muon decay histories
	341
	342	AliMUONTrackLight *mu1 = &fMu0;
	343	AliMUONTrackLight *mu2 = &fMu1;
	344
	345	printf("========= History =======================\n");
	346	printf("first muon");
	347	mu1->PrintInfo("H");
	348	printf("second muon");
	349	mu2->PrintInfo("H");
	350	printf("=========================================\n");
	351	}
	352	if(options.Contains("F") \|\| options.Contains("A")){//flags
	353	printf("the flags set for this muon pair are:\n");
	354	printf("=====================================\n");
	355	if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
	356	fIsCorrelated ? printf("() it is a correlated pair\n") : printf("() it is not a correlated pair\n");
	357	if(IsOpenCharm()) printf("(*) correlated open charm: ");
	358	if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
	359	if(IsOpenCharm() \|\| IsOpenBeauty()){
	360	switch(fCreationProcess){
	361	case 0:
	362	printf("pair creation");
	363	break;
	364	case 1:
	365	printf("gluon splitting");
	366	break;
	367	case 2:
	368	printf("flavour excitation");
	369	break;
	370	case 3:
	371	printf("both muons come from same fragmented mother");
	372	break;
	373	}
	374	if(this->GetMuon(0)->GetOscillation() \|\| this->GetMuon(1)->GetOscillation())
	375	printf("... where oscillation occured\n");
	376	else{
	377	if(IsOpenBeauty())
	378	printf(" (no oscillation)\n");
	379	else
	380	printf("\n");
	381	}
	382	}
b88403f3	383	IsAResonance() ? printf("() it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("() it is not a resonance\n");
b88403f3	384	fIsFeedDown ? printf("() mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(0,fMu0.GetNParents()-2), this->GetMuonMotherPDG(0,fMu0.GetNParents()-1)) : printf("() no feed-down\n");
55fd51b0	385	printf("=====================================\n");
	386	}
	387	if(options.Contains("K") \|\| options.Contains("A")){//dimuon kinematics
	388	Double_t *vtx = this->GetMuon(0)->GetVertex();
	389	TLorentzVector momRec = this->GetPRec();
	390	TLorentzVector momGen = this->GetPGen();
	391	printf("the dimuon charge is %d\n", this->GetCharge());
	392	printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
	393	printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
	394	printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
	395	//rapidity, pT, angles, ...
91e620e2	396	printf("Rec. variables: mass %1.3f, pT %1.3f, pseudo-rapidity %1.3f, openingAngle %1.3f (%1.3f degree), theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
91e620e2	397	momRec.M(), momRec.Pt(), momRec.Eta(),
55fd51b0	398	TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
	399	momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
	400	momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
	401	}
	402	}
	403
b88403f3	404	//====================================
55fd51b0	405	Double_t AliMUONPairLight::GetOpeningAngle() {
b88403f3	406	/// opening angle between the two muons in the lab frame (in degrees)
55fd51b0	407	TLorentzVector pRecMu0 = fMu0.GetPRec();
	408	TLorentzVector pRecMu1 = fMu1.GetPRec();
	409	TVector3 pRecMu03 = pRecMu0.Vect();
	410	TVector3 pRecMu13 = pRecMu1.Vect();
	411	Double_t scalar = pRecMu03.Dot(pRecMu13);
	412	Double_t modMu0 = pRecMu03.Mag();
	413	Double_t modMu1 = pRecMu13.Mag();
	414	Double_t theta = (TMath::ACos(scalar/(modMu0modMu1)))(180./TMath::Pi());
	415	return theta;
	416	}
91e620e2	417	//================================================
	418	Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
	419	///check if we have a correlated pi/K
	420
	421	AliMUONTrackLight mu0 = this->GetMuon(0), mu1 = this->GetMuon(1);
	422	Bool_t fromSameLine = kFALSE;
	423	if (mu0->IsParentPionOrKaon() &&
	424	mu1->IsParentPionOrKaon() &&
	425	mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
	426	) fromSameLine = kTRUE;
	427
	428	return fromSameLine;
	429	}