In AliMUONPairLight, AliMUONTrackLight:
[u/mrichter/AliRoot.git] / MUON / AliMUONPairLight.cxx
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),
61     fMu1(dimuCopy.fMu1), 
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 AliMUONPairLight& AliMUONPairLight::operator=(const AliMUONPairLight& dimuCopy)
86 {
87   // check assignment to self
88   if (this == &dimuCopy) return *this;
89
90   // base class assignment
91   TObject::operator=(dimuCopy);
92
93   // assignment operator
94   fMu0 = dimuCopy.fMu0;
95   fMu1 = dimuCopy.fMu1; 
96   fCreationProcess = dimuCopy.fCreationProcess;
97   fIsCorrelated = dimuCopy.fIsCorrelated; 
98   fCauseOfCorrelation  = dimuCopy.fCauseOfCorrelation;
99   fIsFeedDown = dimuCopy.fIsFeedDown;
100 }
101
102 //====================================
103
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(); 
111
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);
117       return kTRUE;
118     }
119   }
120   return kFALSE; 
121 }
122
123 //====================================
124
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; 
131 }
132
133 //====================================
134
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; } 
140 }
141
142 //====================================
143 void AliMUONPairLight::SetProcess(){
144   /// finds the process related to the muon pair (open charm/beauty, resonance, 
145   /// uncorrelated...) 
146
147   AliMUONTrackLight *mu1 = &fMu0;
148   AliMUONTrackLight *mu2 = &fMu1;
149
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;
166         return;
167       }
168     }
169   }
170
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;
176   }
177
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));
195
196           fCreationProcess = 1; 
197           return;
198         }
199       }
200     }
201   }
202
203   Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
204   Int_t line2 = mu2->GetQuarkPythiaLine(2); 
205
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;
216     }
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;
223       }
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; 
230         return;
231       }
232     }
233   }
234
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;
245     }
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;
252       }
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;
263           return;
264         }
265       }
266     }
267   } 
268
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)){
275
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){
280
281         this->SetCorrelated(kTRUE);
282         fCreationProcess = 1;
283         if(GetCauseOfCorrelation() == -1){
284           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
285         }
286         return;
287       }
288     }
289   }
290
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)){
299
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){
304         
305         this->SetCorrelated(kTRUE);
306         fCreationProcess = 1;
307         if(GetCauseOfCorrelation() == -1){
308           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
309         }
310         return;
311       }
312     }
313   }
314
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)){
323
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){
328         
329         this->SetCorrelated(kTRUE);
330         fCreationProcess = 1;
331         if(GetCauseOfCorrelation() == -1){
332           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
333         }
334         return;
335       }
336     }
337   }
338 }
339
340 //====================================
341 void AliMUONPairLight::SetMuons(const AliMUONTrackLight& mu0, const AliMUONTrackLight& mu1){
342   /// set the two muons 
343   fMu0 = mu0; 
344   fMu1 = mu1; 
345   this->SetProcess();
346
347
348 //====================================
349 void AliMUONPairLight::PrintInfo(const Option_t* opt){
350   /// print information about muon pairs
351   /// Options: 
352   /// - "H" single muons' decay histories
353   /// - "K" dimuon kinematics
354   /// - "F" dimuon flags
355   /// - "A" all variables
356   TString options(opt);
357   options.ToUpper();
358
359   if(options.Contains("H") || options.Contains("A")){//muon decay histories
360
361     AliMUONTrackLight *mu1 = &fMu0;
362     AliMUONTrackLight *mu2 = &fMu1;
363
364     printf("========= History =======================\n");
365     printf("first muon");
366     mu1->PrintInfo("H");
367     printf("second muon");
368     mu2->PrintInfo("H");
369     printf("=========================================\n");
370   }
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){
380       case 0:
381         printf("pair creation");
382         break;
383       case 1:
384         printf("gluon splitting");
385         break;
386       case 2:
387         printf("flavour excitation");
388         break;
389       case 3:
390         printf("both muons come from same fragmented mother");
391         break;
392       }
393       if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation()) 
394         printf("... where oscillation occured\n");
395       else{
396         if(IsOpenBeauty())
397           printf(" (no oscillation)\n");
398         else
399           printf("\n");
400       }
401     }
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");
405   }
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());
420   }
421 }
422
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());
434   return theta; 
435 }
436 //================================================
437 Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
438   ///check if we have a correlated pi/K
439
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;
446
447   return fromSameLine;
448 }