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. *
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 **************************************************************************/
18 //-----------------------------------------------------------------------------
19 //This class was prepared by INFN Cagliari, July 2006
20 //(authors: H.Woehri, A.de Falco)
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
31 // To be used together with AliMUONTrackLight
32 //-----------------------------------------------------------------------------
36 #include "AliMUONPairLight.h"
40 ClassImp(AliMUONPairLight)
42 //====================================
43 AliMUONPairLight::AliMUONPairLight() :
47 fCreationProcess(-999),
48 fIsCorrelated(kFALSE),
49 fCauseOfCorrelation (-1),
52 /// default constructor
56 //====================================
58 AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy)
62 fCreationProcess(dimuCopy.fCreationProcess),
63 fIsCorrelated(dimuCopy.fIsCorrelated),
64 fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
65 fIsFeedDown(dimuCopy.fIsFeedDown)
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;
77 //====================================
79 AliMUONPairLight::~AliMUONPairLight(){
83 //====================================
85 AliMUONPairLight& AliMUONPairLight::operator=(const AliMUONPairLight& dimuCopy)
87 // check assignment to self
88 if (this == &dimuCopy) return *this;
90 // base class assignment
91 TObject::operator=(dimuCopy);
93 // assignment operator
96 fCreationProcess = dimuCopy.fCreationProcess;
97 fIsCorrelated = dimuCopy.fIsCorrelated;
98 fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
99 fIsFeedDown = dimuCopy.fIsFeedDown;
102 //====================================
104 Bool_t AliMUONPairLight::IsAResonance(){
105 /// checks if muon pair comes from a resonance decay
106 if (!fIsCorrelated) return kFALSE; //if muons not correlated, cannot be a resonance
107 //if muons are correlated, check if the PDG of the
108 //common mother is a resonance
109 Int_t nparents0 = fMu0.GetNParents();
110 Int_t nparents1 = fMu1.GetNParents();
112 Int_t minP = TMath::Min(nparents0, nparents1);
113 for (Int_t i = 0 ; i < minP; i++) {
114 if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) &&
115 fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
116 if (nparents0-1-i) SetFeedDown(nparents0-1-i);
123 //====================================
125 AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index) {
126 /// return muon 0 or 1
127 if (index==0) return &fMu0;
128 else if (index==1) return &fMu1;
129 else{ printf ("Index can be either 0 or 1\n"); return 0;}
130 // else return &fMu1;
133 //====================================
135 Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
136 /// return muon mother pdg code
137 if (imuon==0) return fMu0.GetParentPDGCode(mother);
138 else if (imuon==1) return fMu1.GetParentPDGCode(mother);
139 else { printf ("Index must be only 0 or 1\n"); return -999; }
142 //====================================
143 void AliMUONPairLight::SetProcess(){
144 /// finds the process related to the muon pair (open charm/beauty, resonance,
147 AliMUONTrackLight *mu1 = &fMu0;
148 AliMUONTrackLight *mu2 = &fMu1;
150 // check if the two muons are correlated
151 // first check if they come from the same hadron (resonance or beauty/charm meson)
152 Int_t npar1 = mu1->GetNParents();
153 Int_t npar2 = mu2->GetNParents();
154 for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
155 Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
156 for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
157 Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
158 if(lineMo1 == lineMo2) {
159 //reject "diquark" mothers
160 if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
161 // if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
162 this->SetCorrelated(kTRUE);
163 this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
164 if(!IsAResonance()) fCreationProcess = 3;
165 else fCreationProcess = -1;
171 //now, check if we have a correlated pi/K:
172 if(this->IsDimuonFromCorrPiK()){
173 this->SetCorrelated(kTRUE);
174 this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
175 fCreationProcess = -1;
178 // if Open Beauty/Charm we can have 3 creation processes
179 // (pair creation [0], gluon splitting [1] or flavour excitation [2])
180 // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
181 Int_t flavPar1 = mu1->GetParentFlavour(0);
182 Int_t flavPar2 = mu2->GetParentFlavour(0);
183 for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) {
184 Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
185 for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) {
186 Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
187 if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
188 //now, check also that the string fragmented into two hadrons
189 //of the same flavour (string usually splits into many hadrons
190 //among which there are mostly soft particles)
191 if(flavPar1 == flavPar2){
192 this->SetCorrelated(kTRUE);
193 if(GetCauseOfCorrelation() == -1)
194 this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
196 fCreationProcess = 1;
203 Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
204 Int_t line2 = mu2->GetQuarkPythiaLine(2);
206 Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
207 Int_t flavourLine6or7[2] = {-1, -1};
208 // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
209 // are filled with a Q and Qbar
210 for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
211 Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
212 Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
213 if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
214 line6or7[0] = imoth1;
215 flavourLine6or7[0] = flavour1;
217 for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
218 Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
219 Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
220 if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
221 line6or7[1] = imoth2;
222 flavourLine6or7[1] = flavour2;
224 if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
225 (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
226 (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
227 (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour
228 this->SetCorrelated(kTRUE);
229 fCreationProcess = 0;
235 // 3.)flavour excitation: if pythia line 6 of one track *and* pythia line 7 of second track
236 // are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
237 Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
238 Int_t flavourLine2or3[2] = {-1, -1};
239 for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
240 Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
241 Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
242 if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
243 line2or3[0] = imoth1;
244 flavourLine2or3[0] = flavour1;
246 for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
247 Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
248 Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
249 if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
250 line2or3[1] = imoth2;
251 flavourLine2or3[1] = flavour2;
253 if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
254 (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
255 ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) && //or the same,
256 (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
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 this->SetCorrelated(kTRUE);
262 fCreationProcess = 2;
269 //now flag (rare) processes in which only the incoming parton in line 2 or 3
270 //radiates a gluon which produces a QQbar pair:
271 //exclude the light quarks
272 if(line1 == line2 && (line1 == 2 || line1 == 3)){
273 if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
274 (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
276 //now, check also that the string fragmented into two hadrons
277 //of the same flavour (string usually splits into many hadrons
278 //among which there are mostly soft particles)
279 if(flavPar1 == flavPar2){
281 this->SetCorrelated(kTRUE);
282 fCreationProcess = 1;
283 if(GetCauseOfCorrelation() == -1){
284 this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
291 //in initial-state-radiation produced QQbar events the "mother quark"
292 //is acknowledged as the second quark [1] and sits in line 2 or 3
293 //is part of gluon splitting
294 line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
295 line2 = mu2->GetQuarkPythiaLine(1);
296 if(line1 == line2 && (line1 == 2 || line1 == 3)){
297 if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
298 (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
300 //now, check also that the string fragmented into two hadrons
301 //of the same flavour (string usually splits into many hadrons
302 //among which there are mostly soft particles)
303 if(flavPar1 == flavPar2){
305 this->SetCorrelated(kTRUE);
306 fCreationProcess = 1;
307 if(GetCauseOfCorrelation() == -1){
308 this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
315 //in final-state-radiation produced QQbar events the "mother quark"
316 //is acknowledged as the first quark [1] and sits in line 6 or 7
317 //is part of gluon splitting
318 line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark
319 line2 = mu2->GetQuarkPythiaLine(1);
320 if(line1 == line2 && (line1 == 6 || line1 == 7)){
321 if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
322 (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
324 //now, check also that the string fragmented into two hadrons
325 //of the same flavour (string usually splits into many hadrons
326 //among which there are mostly soft particles)
327 if(flavPar1 == flavPar2){
329 this->SetCorrelated(kTRUE);
330 fCreationProcess = 1;
331 if(GetCauseOfCorrelation() == -1){
332 this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
340 //====================================
341 void AliMUONPairLight::SetMuons(const AliMUONTrackLight& mu0, const AliMUONTrackLight& mu1){
342 /// set the two muons
348 //====================================
349 void AliMUONPairLight::PrintInfo(const Option_t* opt){
350 /// print information about muon pairs
352 /// - "H" single muons' decay histories
353 /// - "K" dimuon kinematics
354 /// - "F" dimuon flags
355 /// - "A" all variables
356 TString options(opt);
359 if(options.Contains("H") || options.Contains("A")){//muon decay histories
361 AliMUONTrackLight *mu1 = &fMu0;
362 AliMUONTrackLight *mu2 = &fMu1;
364 printf("========= History =======================\n");
365 printf("first muon");
367 printf("second muon");
369 printf("=========================================\n");
371 if(options.Contains("F") || options.Contains("A")){//flags
372 printf("the flags set for this muon pair are:\n");
373 printf("=====================================\n");
374 if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
375 fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
376 if(IsOpenCharm()) printf("(*) correlated open charm: ");
377 if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
378 if(IsOpenCharm() || IsOpenBeauty()){
379 switch(fCreationProcess){
381 printf("pair creation");
384 printf("gluon splitting");
387 printf("flavour excitation");
390 printf("both muons come from same fragmented mother");
393 if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation())
394 printf("... where oscillation occured\n");
397 printf(" (no oscillation)\n");
402 IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
403 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");
404 printf("=====================================\n");
406 if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
407 Double_t *vtx = this->GetMuon(0)->GetVertex();
408 TLorentzVector momRec = this->GetPRec();
409 TLorentzVector momGen = this->GetPGen();
410 printf("the dimuon charge is %d\n", this->GetCharge());
411 printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
412 printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
413 printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
414 //rapidity, pT, angles, ...
415 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",
416 momRec.M(), momRec.Pt(), momRec.Eta(),
417 TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
418 momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
419 momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
423 //====================================
424 Double_t AliMUONPairLight::GetOpeningAngle() {
425 /// opening angle between the two muons in the lab frame (in degrees)
426 TLorentzVector pRecMu0 = fMu0.GetPRec();
427 TLorentzVector pRecMu1 = fMu1.GetPRec();
428 TVector3 pRecMu03 = pRecMu0.Vect();
429 TVector3 pRecMu13 = pRecMu1.Vect();
430 Double_t scalar = pRecMu03.Dot(pRecMu13);
431 Double_t modMu0 = pRecMu03.Mag();
432 Double_t modMu1 = pRecMu13.Mag();
433 Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
436 //================================================
437 Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
438 ///check if we have a correlated pi/K
440 AliMUONTrackLight *mu0 = this->GetMuon(0), *mu1 = this->GetMuon(1);
441 Bool_t fromSameLine = kFALSE;
442 if (mu0->IsParentPionOrKaon() &&
443 mu1->IsParentPionOrKaon() &&
444 mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
445 ) fromSameLine = kTRUE;