bug fixed
[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
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
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),
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
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
b88403f3 90 Int_t nparents0 = fMu0.GetNParents();
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
106AliMUONTrackLight* 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
116Int_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//====================================
124void 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--) {
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 }
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
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
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
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//====================================
322void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
323 /// set the two muons
324 fMu0 = mu0;
325 fMu1 = mu1;
326 this->SetProcess();
327}
328
329//====================================
57e2ad1a 330void AliMUONPairLight::PrintInfo(const Option_t* opt){
55fd51b0 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");
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",
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 405Double_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/(modMu0*modMu1)))*(180./TMath::Pi());
415 return theta;
416}
91e620e2 417//================================================
418Bool_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}