Line data Source code
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 16 : ClassImp(AliMUONPairLight)
41 :
42 : //====================================
43 : AliMUONPairLight::AliMUONPairLight() :
44 0 : TObject(),
45 0 : fMu0(),
46 0 : fMu1(),
47 0 : fCreationProcess(-999),
48 0 : fIsCorrelated(kFALSE),
49 0 : fCauseOfCorrelation (-1),
50 0 : fIsFeedDown(kFALSE)
51 0 : {
52 : /// default constructor
53 : ;
54 0 : }
55 :
56 : //====================================
57 :
58 : AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy)
59 0 : : TObject(dimuCopy),
60 0 : fMu0(dimuCopy.fMu0),
61 0 : fMu1(dimuCopy.fMu1),
62 0 : fCreationProcess(dimuCopy.fCreationProcess),
63 0 : fIsCorrelated(dimuCopy.fIsCorrelated),
64 0 : fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
65 0 : fIsFeedDown(dimuCopy.fIsFeedDown)
66 0 : {
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 0 : }
76 :
77 : //====================================
78 :
79 0 : AliMUONPairLight::~AliMUONPairLight(){
80 : /// destructor
81 0 : }
82 :
83 : //====================================
84 :
85 : AliMUONPairLight& AliMUONPairLight::operator=(const AliMUONPairLight& dimuCopy)
86 : {
87 : // check assignment to self
88 0 : if (this == &dimuCopy) return *this;
89 :
90 : // base class assignment
91 0 : TObject::operator=(dimuCopy);
92 :
93 : // assignment operator
94 0 : fMu0 = dimuCopy.fMu0;
95 0 : fMu1 = dimuCopy.fMu1;
96 0 : fCreationProcess = dimuCopy.fCreationProcess;
97 0 : fIsCorrelated = dimuCopy.fIsCorrelated;
98 0 : fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
99 0 : fIsFeedDown = dimuCopy.fIsFeedDown;
100 :
101 0 : return *this;
102 0 : }
103 :
104 : //====================================
105 :
106 : Bool_t AliMUONPairLight::IsAResonance(){
107 : /// checks if muon pair comes from a resonance decay
108 0 : if (!fIsCorrelated) return kFALSE; //if muons not correlated, cannot be a resonance
109 : //if muons are correlated, check if the PDG of the
110 : //common mother is a resonance
111 0 : Int_t nparents0 = fMu0.GetNParents();
112 0 : Int_t nparents1 = fMu1.GetNParents();
113 :
114 0 : Int_t minP = TMath::Min(nparents0, nparents1);
115 0 : for (Int_t i = 0 ; i < minP; i++) {
116 0 : if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) &&
117 0 : fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
118 0 : if (nparents0-1-i) SetFeedDown(nparents0-1-i);
119 0 : return kTRUE;
120 : }
121 : }
122 0 : return kFALSE;
123 0 : }
124 :
125 : //====================================
126 :
127 : AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index) {
128 : /// return muon 0 or 1
129 0 : if (index==0) return &fMu0;
130 0 : else if (index==1) return &fMu1;
131 0 : else{ printf ("Index can be either 0 or 1\n"); return 0;}
132 : // else return &fMu1;
133 0 : }
134 :
135 : //====================================
136 :
137 : Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
138 : /// return muon mother pdg code
139 0 : if (imuon==0) return fMu0.GetParentPDGCode(mother);
140 0 : else if (imuon==1) return fMu1.GetParentPDGCode(mother);
141 0 : else { printf ("Index must be only 0 or 1\n"); return -999; }
142 0 : }
143 :
144 : //====================================
145 : void AliMUONPairLight::SetProcess(){
146 : /// finds the process related to the muon pair (open charm/beauty, resonance,
147 : /// uncorrelated...)
148 :
149 0 : AliMUONTrackLight *mu1 = &fMu0;
150 0 : AliMUONTrackLight *mu2 = &fMu1;
151 :
152 : // check if the two muons are correlated
153 : // first check if they come from the same hadron (resonance or beauty/charm meson)
154 0 : Int_t npar1 = mu1->GetNParents();
155 0 : Int_t npar2 = mu2->GetNParents();
156 0 : for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
157 0 : Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
158 0 : for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
159 0 : Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
160 0 : if(lineMo1 == lineMo2) {
161 : //reject "diquark" mothers
162 0 : if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
163 : // if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
164 0 : this->SetCorrelated(kTRUE);
165 0 : this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
166 0 : if(!IsAResonance()) fCreationProcess = 3;
167 0 : else fCreationProcess = -1;
168 0 : return;
169 : }
170 0 : }
171 0 : }
172 :
173 : //now, check if we have a correlated pi/K:
174 0 : if(this->IsDimuonFromCorrPiK()){
175 0 : this->SetCorrelated(kTRUE);
176 0 : this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
177 0 : fCreationProcess = -1;
178 0 : }
179 :
180 : // if Open Beauty/Charm we can have 3 creation processes
181 : // (pair creation [0], gluon splitting [1] or flavour excitation [2])
182 : // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
183 0 : Int_t flavPar1 = mu1->GetParentFlavour(0);
184 0 : Int_t flavPar2 = mu2->GetParentFlavour(0);
185 0 : for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) {
186 0 : Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
187 0 : for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) {
188 0 : Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
189 0 : if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
190 : //now, check also that the string fragmented into two hadrons
191 : //of the same flavour (string usually splits into many hadrons
192 : //among which there are mostly soft particles)
193 0 : if(flavPar1 == flavPar2){
194 0 : this->SetCorrelated(kTRUE);
195 0 : if(GetCauseOfCorrelation() == -1)
196 0 : this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
197 :
198 0 : fCreationProcess = 1;
199 0 : return;
200 : }
201 : }
202 0 : }
203 0 : }
204 :
205 0 : Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
206 0 : Int_t line2 = mu2->GetQuarkPythiaLine(2);
207 :
208 : Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
209 : Int_t flavourLine6or7[2] = {-1, -1};
210 : // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
211 : // are filled with a Q and Qbar
212 0 : for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
213 0 : Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
214 0 : Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
215 0 : if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
216 : line6or7[0] = imoth1;
217 : flavourLine6or7[0] = flavour1;
218 0 : }
219 0 : for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
220 0 : Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
221 0 : Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
222 0 : if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
223 : line6or7[1] = imoth2;
224 : flavourLine6or7[1] = flavour2;
225 0 : }
226 0 : if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
227 0 : (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
228 0 : (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
229 0 : (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour
230 0 : this->SetCorrelated(kTRUE);
231 0 : fCreationProcess = 0;
232 0 : return;
233 : }
234 0 : }
235 0 : }
236 :
237 : // 3.)flavour excitation: if pythia line 6 of one track *and* pythia line 7 of second track
238 : // are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
239 : Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
240 : Int_t flavourLine2or3[2] = {-1, -1};
241 0 : for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
242 0 : Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
243 0 : Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
244 0 : if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
245 : line2or3[0] = imoth1;
246 : flavourLine2or3[0] = flavour1;
247 0 : }
248 0 : for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
249 0 : Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
250 0 : Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
251 0 : if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
252 : line2or3[1] = imoth2;
253 : flavourLine2or3[1] = flavour2;
254 0 : }
255 0 : if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
256 0 : (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
257 0 : ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) && //or the same,
258 0 : (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
259 : //now, check also that the string fragmented into two hadrons
260 : //of the same flavour (string usually splits into many hadrons
261 : //among which there are mostly soft particles)
262 0 : if(flavPar1 == flavPar2){
263 0 : this->SetCorrelated(kTRUE);
264 0 : fCreationProcess = 2;
265 0 : return;
266 : }
267 : }
268 0 : }
269 0 : }
270 :
271 : //now flag (rare) processes in which only the incoming parton in line 2 or 3
272 : //radiates a gluon which produces a QQbar pair:
273 : //exclude the light quarks
274 0 : if(line1 == line2 && (line1 == 2 || line1 == 3)){
275 0 : if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
276 0 : (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
277 :
278 : //now, check also that the string fragmented into two hadrons
279 : //of the same flavour (string usually splits into many hadrons
280 : //among which there are mostly soft particles)
281 0 : if(flavPar1 == flavPar2){
282 :
283 0 : this->SetCorrelated(kTRUE);
284 0 : fCreationProcess = 1;
285 0 : if(GetCauseOfCorrelation() == -1){
286 0 : this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
287 0 : }
288 0 : return;
289 : }
290 : }
291 : }
292 :
293 : //in initial-state-radiation produced QQbar events the "mother quark"
294 : //is acknowledged as the second quark [1] and sits in line 2 or 3
295 : //is part of gluon splitting
296 0 : line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
297 0 : line2 = mu2->GetQuarkPythiaLine(1);
298 0 : if(line1 == line2 && (line1 == 2 || line1 == 3)){
299 0 : if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
300 0 : (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
301 :
302 : //now, check also that the string fragmented into two hadrons
303 : //of the same flavour (string usually splits into many hadrons
304 : //among which there are mostly soft particles)
305 0 : if(flavPar1 == flavPar2){
306 :
307 0 : this->SetCorrelated(kTRUE);
308 0 : fCreationProcess = 1;
309 0 : if(GetCauseOfCorrelation() == -1){
310 0 : this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
311 0 : }
312 0 : return;
313 : }
314 : }
315 : }
316 :
317 : //in final-state-radiation produced QQbar events the "mother quark"
318 : //is acknowledged as the first quark [1] and sits in line 6 or 7
319 : //is part of gluon splitting
320 0 : line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark
321 0 : line2 = mu2->GetQuarkPythiaLine(1);
322 0 : if(line1 == line2 && (line1 == 6 || line1 == 7)){
323 0 : if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
324 0 : (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
325 :
326 : //now, check also that the string fragmented into two hadrons
327 : //of the same flavour (string usually splits into many hadrons
328 : //among which there are mostly soft particles)
329 0 : if(flavPar1 == flavPar2){
330 :
331 0 : this->SetCorrelated(kTRUE);
332 0 : fCreationProcess = 1;
333 0 : if(GetCauseOfCorrelation() == -1){
334 0 : this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
335 0 : }
336 0 : return;
337 : }
338 : }
339 : }
340 0 : }
341 :
342 : //====================================
343 : void AliMUONPairLight::SetMuons(const AliMUONTrackLight& mu0, const AliMUONTrackLight& mu1){
344 : /// set the two muons
345 0 : fMu0 = mu0;
346 0 : fMu1 = mu1;
347 0 : this->SetProcess();
348 0 : }
349 :
350 : //====================================
351 : void AliMUONPairLight::PrintInfo(const Option_t* opt){
352 : /// print information about muon pairs
353 : /// Options:
354 : /// - "H" single muons' decay histories
355 : /// - "K" dimuon kinematics
356 : /// - "F" dimuon flags
357 : /// - "A" all variables
358 0 : TString options(opt);
359 0 : options.ToUpper();
360 :
361 0 : if(options.Contains("H") || options.Contains("A")){//muon decay histories
362 :
363 0 : AliMUONTrackLight *mu1 = &fMu0;
364 0 : AliMUONTrackLight *mu2 = &fMu1;
365 :
366 0 : printf("========= History =======================\n");
367 0 : printf("first muon");
368 0 : mu1->PrintInfo("H");
369 0 : printf("second muon");
370 0 : mu2->PrintInfo("H");
371 0 : printf("=========================================\n");
372 0 : }
373 0 : if(options.Contains("F") || options.Contains("A")){//flags
374 0 : printf("the flags set for this muon pair are:\n");
375 0 : printf("=====================================\n");
376 0 : if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
377 0 : fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
378 0 : if(IsOpenCharm()) printf("(*) correlated open charm: ");
379 0 : if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
380 0 : if(IsOpenCharm() || IsOpenBeauty()){
381 0 : switch(fCreationProcess){
382 : case 0:
383 0 : printf("pair creation");
384 : break;
385 : case 1:
386 0 : printf("gluon splitting");
387 : break;
388 : case 2:
389 0 : printf("flavour excitation");
390 : break;
391 : case 3:
392 0 : printf("both muons come from same fragmented mother");
393 : break;
394 : }
395 0 : if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation())
396 0 : printf("... where oscillation occured\n");
397 : else{
398 0 : if(IsOpenBeauty())
399 0 : printf(" (no oscillation)\n");
400 : else
401 0 : printf("\n");
402 : }
403 : }
404 0 : IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
405 0 : 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");
406 0 : printf("=====================================\n");
407 : }
408 0 : if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
409 0 : Double_t *vtx = this->GetMuon(0)->GetVertex();
410 0 : TLorentzVector momRec = this->GetPRec();
411 0 : TLorentzVector momGen = this->GetPGen();
412 0 : printf("the dimuon charge is %d\n", this->GetCharge());
413 0 : printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
414 0 : printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
415 0 : printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
416 : //rapidity, pT, angles, ...
417 0 : 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",
418 0 : momRec.M(), momRec.Pt(), momRec.Eta(),
419 0 : TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
420 0 : momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
421 0 : momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
422 0 : }
423 0 : }
424 :
425 : //====================================
426 : Double_t AliMUONPairLight::GetOpeningAngle() {
427 : /// opening angle between the two muons in the lab frame (in degrees)
428 0 : TLorentzVector pRecMu0 = fMu0.GetPRec();
429 0 : TLorentzVector pRecMu1 = fMu1.GetPRec();
430 0 : TVector3 pRecMu03 = pRecMu0.Vect();
431 0 : TVector3 pRecMu13 = pRecMu1.Vect();
432 0 : Double_t scalar = pRecMu03.Dot(pRecMu13);
433 0 : Double_t modMu0 = pRecMu03.Mag();
434 0 : Double_t modMu1 = pRecMu13.Mag();
435 0 : Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
436 : return theta;
437 0 : }
438 : //================================================
439 : Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
440 : ///check if we have a correlated pi/K
441 :
442 0 : AliMUONTrackLight *mu0 = this->GetMuon(0), *mu1 = this->GetMuon(1);
443 : Bool_t fromSameLine = kFALSE;
444 0 : if (mu0->IsParentPionOrKaon() &&
445 0 : mu1->IsParentPionOrKaon() &&
446 0 : mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
447 0 : ) fromSameLine = kTRUE;
448 :
449 0 : return fromSameLine;
450 : }
|
