[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 pdg = GetCauseOfCorrelation();
  if(pdg < 10) return kFALSE;
  Int_t id=pdg%100000;
  if(((id-id%10)%110)) return kFALSE;
  else return kTRUE;
  printf("<AliMUONPairLight::IsAResonance> arriving after this piece of code\n");

}

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

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;

  //1.) check if one of the tracks is not a muon
  if(IsOneTrackNotAMuon()) { 
    this->SetCorrelated(kFALSE);
    return;
  }

  // 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 = 0; imoth1<npar1; imoth1++) { 
  for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) { 
    Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
    //    for (Int_t imoth2 = 0; imoth2<npar2; imoth2++) { 
    for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) { 
      Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
      if(lineMo1 == lineMo2) { 
	this->SetCorrelated(kTRUE); 
	this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
	if(!IsAResonance()) fCreationProcess = 3; 
	else fCreationProcess = -1;
	return;   // <<<<---------------- RETURN? 
      }
    }
  }  
  
  // 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 
  if (mu1->GetQuarkPythiaLine(2) == mu2->GetQuarkPythiaLine(2) && mu1->GetQuarkPDGCode(2) == 21) { 
    this->SetCorrelated(kTRUE); 
    if(GetCauseOfCorrelation() == -1){
      this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(2));
    }
    fCreationProcess = 1; 
    return;
  }
  
  // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
  // are filled with a Q and Qbar
  Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
  Int_t line2 = mu2->GetQuarkPythiaLine(2); 
  Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(2)); //[2] ... very first quark
  Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(2)); 
  if ((line1 == 6 || line1 == 7) && (line2 == 6 || line2 == 7)) { 
    if((flavour1 == 4 || flavour1 == 5) && (flavour2 == 4 || flavour2 == 5)){
      this->SetCorrelated(kTRUE);
      fCreationProcess = 0; 
      return;
    }
  }

  // 3.)flavour excitation:
  if((((line1 == 6 || line1 == 7) && (flavour1 == 4 || flavour1 == 5)) && //first track has Q in line 6 or 7
      ((line2 == 2 || line2 == 3) && (flavour2 == 21 || flavour2 < 10))) || //second track has a g/q in line 2 or 3
      (((line2 == 6 || line2 == 7) && (flavour2 == 4 || flavour2 == 5)) && //or the same,
      ((line1 == 2 || line1 == 3) && (flavour1 == 21 || flavour1 < 10)))){ // swapping the track's indices

//     printf("candidate for flavour excitation\n");

    //Hermine: is it needed to check the lines 4-7???

    //now, we have a candidate for flavour excitation
    //we must also verify if in the Pythia listing
    //the "incoming" lines 4 and 5 and "outgoing" lines 6 and 7
    //are filled with g Q each: e.g. 4(g),5(Q),5(g),6(Q)
//     Int_t pdg4 = TMath::Abs(stack()->Particle(4)->GetPdgCode());
//     Int_t pdg5 = TMath::Abs(stack()->Particle(5)->GetPdgCode());
//     Int_t pdg6 = TMath::Abs(stack()->Particle(6)->GetPdgCode());
//     Int_t pdg7 = TMath::Abs(stack()->Particle(7)->GetPdgCode());
//     if((pdg4 == 21 && pdg5 < 10 || pdg5 == 21 && pdg4 < 10 ) &&
//        (pdg6 == 21 && pdg7 < 10 || pdg7 == 21 && pdg6 < 10 )){
//       this->PrintInfo("H");
      this->SetCorrelated(kTRUE);
      fCreationProcess = 2;
      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(fIsFeedDown)) : printf("(*) it is not a resonance\n");
    fIsFeedDown ? printf("(*) mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(1), this->GetMuonMotherPDG(0)) : 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: 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.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; 
}
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
	18	//===================================================================
	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
	32	//===================================================================
	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
	90	Int_t pdg = GetCauseOfCorrelation();
	91	if(pdg < 10) return kFALSE;
	92	Int_t id=pdg%100000;
	93	if(((id-id%10)%110)) return kFALSE;
	94	else return kTRUE;
	95	printf("<AliMUONPairLight::IsAResonance> arriving after this piece of code\n");
	96
	97	}
	98
	99	//====================================
	100
	101	AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index) {
	102	/// return muon 0 or 1
	103	if (index==0) return &fMu0;
	104	else if (index==1) return &fMu1;
	105	else{ printf ("Index can be either 0 or 1\n"); return 0;}
	106	// else return &fMu1;
	107	}
	108
	109	//====================================
	110
	111	Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
	112	/// return muon mother pdg code
	113	if (imuon==0) return fMu0.GetParentPDGCode(mother);
	114	else if (imuon==1) return fMu1.GetParentPDGCode(mother);
	115	else { printf ("Index must be only 0 or 1\n"); return -999; }
	116	}
	117
	118	//====================================
	119	void AliMUONPairLight::SetProcess(){
	120	/// finds the process related to the muon pair (open charm/beauty, resonance,
	121	/// uncorrelated...)
	122	AliMUONTrackLight *mu1 = &fMu0;
	123	AliMUONTrackLight *mu2 = &fMu1;
	124
	125	//1.) check if one of the tracks is not a muon
126	if(IsOneTrackNotAMuon()) {
127	this->SetCorrelated(kFALSE);
128	return;
129	}
130
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();
135	// for (Int_t imoth1 = 0; imoth1<npar1; imoth1++) {
136	for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
137	Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
138	// for (Int_t imoth2 = 0; imoth2<npar2; imoth2++) {
139	for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
140	Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
141	if(lineMo1 == lineMo2) {
142	this->SetCorrelated(kTRUE);
143	this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
144	if(!IsAResonance()) fCreationProcess = 3;
145	else fCreationProcess = -1;
146	return; // <<<<---------------- RETURN?
147	}
148	}
149	}
150
151	// if Open Beauty/Charm we can have 3 creation processes
152	// (pair creation [0], gluon splitting [1] or flavour excitation [2])
153	//
154	// 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
155	if (mu1->GetQuarkPythiaLine(2) == mu2->GetQuarkPythiaLine(2) && mu1->GetQuarkPDGCode(2) == 21) {
156	this->SetCorrelated(kTRUE);
157	if(GetCauseOfCorrelation() == -1){
158	this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(2));
159	}
160	fCreationProcess = 1;
161	return;
162	}
163
164	// 2.) pair creation: if pythia line 6 of one track and pythia line 7 of second track
165	// are filled with a Q and Qbar
166	Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
167	Int_t line2 = mu2->GetQuarkPythiaLine(2);
168	Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(2)); //[2] ... very first quark
169	Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(2));
170	if ((line1 == 6 \|\| line1 == 7) && (line2 == 6 \|\| line2 == 7)) {
171	if((flavour1 == 4 \|\| flavour1 == 5) && (flavour2 == 4 \|\| flavour2 == 5)){
172	this->SetCorrelated(kTRUE);
173	fCreationProcess = 0;
174	return;
175	}
176	}
177
178	// 3.)flavour excitation:
179	if((((line1 == 6 \|\| line1 == 7) && (flavour1 == 4 \|\| flavour1 == 5)) && //first track has Q in line 6 or 7
180	((line2 == 2 \|\| line2 == 3) && (flavour2 == 21 \|\| flavour2 < 10))) \|\| //second track has a g/q in line 2 or 3
181	(((line2 == 6 \|\| line2 == 7) && (flavour2 == 4 \|\| flavour2 == 5)) && //or the same,
182	((line1 == 2 \|\| line1 == 3) && (flavour1 == 21 \|\| flavour1 < 10)))){ // swapping the track's indices
183
184	// printf("candidate for flavour excitation\n");
185
186	//Hermine: is it needed to check the lines 4-7???
187
188	//now, we have a candidate for flavour excitation
189	//we must also verify if in the Pythia listing
190	//the "incoming" lines 4 and 5 and "outgoing" lines 6 and 7
191	//are filled with g Q each: e.g. 4(g),5(Q),5(g),6(Q)
192	// Int_t pdg4 = TMath::Abs(stack()->Particle(4)->GetPdgCode());
193	// Int_t pdg5 = TMath::Abs(stack()->Particle(5)->GetPdgCode());
194	// Int_t pdg6 = TMath::Abs(stack()->Particle(6)->GetPdgCode());
195	// Int_t pdg7 = TMath::Abs(stack()->Particle(7)->GetPdgCode());
196	// if((pdg4 == 21 && pdg5 < 10 \|\| pdg5 == 21 && pdg4 < 10 ) &&
197	// (pdg6 == 21 && pdg7 < 10 \|\| pdg7 == 21 && pdg6 < 10 )){
198	// this->PrintInfo("H");
199	this->SetCorrelated(kTRUE);
200	fCreationProcess = 2;
201	return;
202	// }
203	}
204	}
205
206	//====================================
207	void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
208	/// set the two muons
209	fMu0 = mu0;
210	fMu1 = mu1;
211	this->SetProcess();
212	}
213
214	//====================================
215	void AliMUONPairLight::PrintInfo(Option_t* opt){
216	/// print information about muon pairs
217	/// Options:
218	/// - "H" single muons' decay histories
219	/// - "K" dimuon kinematics
220	/// - "F" dimuon flags
221	/// - "A" all variables
222	TString options(opt);
223	options.ToUpper();
224
225	if(options.Contains("H") \|\| options.Contains("A")){//muon decay histories
226
227	AliMUONTrackLight *mu1 = &fMu0;
228	AliMUONTrackLight *mu2 = &fMu1;
229
230	printf("========= History =======================\n");
231	printf("first muon");
232	mu1->PrintInfo("H");
233	printf("second muon");
234	mu2->PrintInfo("H");
235	printf("=========================================\n");
236	}
237	if(options.Contains("F") \|\| options.Contains("A")){//flags
238	printf("the flags set for this muon pair are:\n");
239	printf("=====================================\n");
240	if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
241	fIsCorrelated ? printf("() it is a correlated pair\n") : printf("() it is not a correlated pair\n");
242	if(IsOpenCharm()) printf("(*) correlated open charm: ");
243	if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
244	if(IsOpenCharm() \|\| IsOpenBeauty()){
245	switch(fCreationProcess){
246	case 0:
247	printf("pair creation");
248	break;
249	case 1:
250	printf("gluon splitting");
251	break;
252	case 2:
253	printf("flavour excitation");
254	break;
255	case 3:
256	printf("both muons come from same fragmented mother");
257	break;
258	}
259	if(this->GetMuon(0)->GetOscillation() \|\| this->GetMuon(1)->GetOscillation())
260	printf("... where oscillation occured\n");
261	else{
262	if(IsOpenBeauty())
263	printf(" (no oscillation)\n");
264	else
265	printf("\n");
266	}
267	}
268	IsAResonance() ? printf("() it is a resonance: %d\n", this->GetMuonMotherPDG(fIsFeedDown)) : printf("() it is not a resonance\n");
269	fIsFeedDown ? printf("() mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(1), this->GetMuonMotherPDG(0)) : printf("() no feed-down\n");
270	printf("=====================================\n");
271	}
272	if(options.Contains("K") \|\| options.Contains("A")){//dimuon kinematics
273	Double_t *vtx = this->GetMuon(0)->GetVertex();
274	TLorentzVector momRec = this->GetPRec();
275	TLorentzVector momGen = this->GetPGen();
276	printf("the dimuon charge is %d\n", this->GetCharge());
277	printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
278	printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
279	printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
280	//rapidity, pT, angles, ...
281	printf("Rec. variables: 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",
282	momRec.Pt(), momRec.Eta(),
283	TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
284	momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
285	momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
286	}
287	}
288
289	Double_t AliMUONPairLight::GetOpeningAngle() {
290	// opening angle between the two muons in the lab frame (in degrees)
291	TLorentzVector pRecMu0 = fMu0.GetPRec();
292	TLorentzVector pRecMu1 = fMu1.GetPRec();
293	TVector3 pRecMu03 = pRecMu0.Vect();
294	TVector3 pRecMu13 = pRecMu1.Vect();
295	Double_t scalar = pRecMu03.Dot(pRecMu13);
296	Double_t modMu0 = pRecMu03.Mag();
297	Double_t modMu1 = pRecMu13.Mag();
298	Double_t theta = (TMath::ACos(scalar/(modMu0modMu1)))(180./TMath::Pi());
299	return theta;
300	}