New classes for analysis
[u/mrichter/AliRoot.git] / MUON / AliMUONPairLight.cxx
CommitLineData
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
40ClassImp(AliMUONPairLight)
41
42//====================================
43AliMUONPairLight::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
58AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy)
59 : TObject(dimuCopy),
60 fMu0(dimuCopy.fMu0),
61 fMu1(dimuCopy.fMu0),
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
79AliMUONPairLight::~AliMUONPairLight(){
80 /// destructor
81}
82
83//====================================
84
85Bool_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
101AliMUONTrackLight* 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
111Int_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//====================================
119void 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//====================================
207void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
208 /// set the two muons
209 fMu0 = mu0;
210 fMu1 = mu1;
211 this->SetProcess();
212}
213
214//====================================
215void 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
289Double_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/(modMu0*modMu1)))*(180./TMath::Pi());
299 return theta;
300}