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 : // Compact information for the muon generated tracks in the MUON arm
20 : // useful at the last stage of the analysis chain
21 : // provides a link between the reconstructed track and the generated particle
22 : // stores kinematical information at gen. and rec. level and
23 : // the decay history of the muon, allowing the identification of the
24 : // mother process
25 : //
26 : // To be used together with AliMUONPairLight
27 : //
28 : // This class was prepared by INFN Cagliari, July 2006
29 : // (authors: H.Woehri, A.de Falco)
30 : //-----------------------------------------------------------------------------
31 :
32 : // 13 Nov 2007:
33 : // Added a temporary fix to FindRefTrack to be able to handle reconstructed tracks
34 : // generated from ESD muon track information. The problem is that the ESD data at
35 : // the moment only contains the first hit on chamber 1. Hopefully in the near future
36 : // this will be fixed and all hit information will be available.
37 : // - Artur Szostak <artursz@iafrica.com>
38 :
39 : #include "AliMUONTrackLight.h"
40 : #include "AliMUONTrack.h"
41 : #include "AliMUONConstants.h"
42 : #include "AliMUONVTrackStore.h"
43 : #include "AliMUONTrackParam.h"
44 :
45 : #include "AliESDMuonTrack.h"
46 : #include "AliStack.h"
47 : #include "AliLog.h"
48 :
49 : #include "TDatabasePDG.h"
50 : #include "TParticle.h"
51 : #include "TString.h"
52 :
53 : #include <cstdio>
54 : #include <Riostream.h>
55 :
56 16 : ClassImp(AliMUONTrackLight)
57 :
58 : //===================================================================
59 :
60 : AliMUONTrackLight::AliMUONTrackLight()
61 0 : : TObject(),
62 0 : fPrec(),
63 0 : fIsTriggered(kFALSE),
64 0 : fCharge(-999),
65 0 : fChi2(-1),
66 0 : fCentr(-1),
67 0 : fPgen(),
68 0 : fTrackPythiaLine(-999),
69 0 : fTrackPDGCode(-999),
70 0 : fOscillation(kFALSE),
71 0 : fNParents(0),
72 0 : fWeight(1)
73 0 : {
74 : /// default constructor
75 0 : fPgen.SetPxPyPzE(0.,0.,0.,0.);
76 0 : fPrec.SetPxPyPzE(0.,0.,0.,0.);
77 0 : for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
78 0 : for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
79 0 : fParentPDGCode[npar] = -1;
80 0 : fParentPythiaLine[npar] = -1;
81 : }
82 0 : for (Int_t i = 0; i < 4; i++){
83 0 : fQuarkPDGCode[i] = -1;
84 0 : fQuarkPythiaLine[i] = -1;
85 : }
86 0 : }
87 :
88 : //============================================
89 : AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
90 0 : : TObject(muonCopy),
91 0 : fPrec(muonCopy.fPrec),
92 0 : fIsTriggered(muonCopy.fIsTriggered),
93 0 : fCharge(muonCopy.fCharge),
94 0 : fChi2(muonCopy.fChi2),
95 0 : fCentr(muonCopy.fCentr),
96 0 : fPgen(muonCopy.fPgen),
97 0 : fTrackPythiaLine(muonCopy.fTrackPythiaLine),
98 0 : fTrackPDGCode(muonCopy.fTrackPDGCode),
99 0 : fOscillation(muonCopy.fOscillation),
100 0 : fNParents(muonCopy.fNParents),
101 0 : fWeight(muonCopy.fWeight)
102 0 : {
103 : /// copy constructor
104 0 : for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
105 0 : for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
106 0 : fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
107 0 : fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
108 : }
109 0 : for (Int_t i = 0; i < 4; i++){
110 0 : fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
111 0 : fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
112 : }
113 0 : }
114 :
115 : //============================================
116 : AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
117 0 : : TObject(),
118 0 : fPrec(),
119 0 : fIsTriggered(kFALSE),
120 0 : fCharge(-999),
121 0 : fChi2(-1),
122 0 : fCentr(-1),
123 0 : fPgen(),
124 0 : fTrackPythiaLine(-999),
125 0 : fTrackPDGCode(-999),
126 0 : fOscillation(kFALSE),
127 0 : fNParents(0),
128 0 : fWeight(1)
129 0 : {
130 : /// constructor
131 0 : fPgen.SetPxPyPzE(0.,0.,0.,0.);
132 0 : for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
133 0 : fParentPDGCode[npar] = -1;
134 0 : fParentPythiaLine[npar] = -1;
135 : }
136 0 : for (Int_t i = 0; i < 4; i++){
137 0 : fQuarkPDGCode[i] = -1;
138 0 : fQuarkPythiaLine[i] = -1;
139 : }
140 0 : FillFromESD(muonTrack);
141 0 : }
142 :
143 : //============================================
144 : AliMUONTrackLight::~AliMUONTrackLight()
145 0 : {
146 : /// Destructor
147 0 : }
148 :
149 : //============================================
150 : AliMUONTrackLight& AliMUONTrackLight::operator=(const AliMUONTrackLight& muonCopy)
151 : {
152 : // check assignment to self
153 0 : if (this == &muonCopy) return *this;
154 :
155 : // base class assignment
156 0 : TObject::operator=(muonCopy);
157 :
158 : // assignment operator
159 0 : fPrec = muonCopy.fPrec;
160 0 : fIsTriggered = muonCopy.fIsTriggered;
161 0 : fCharge = muonCopy.fCharge;
162 0 : fChi2 = muonCopy.fChi2;
163 0 : fCentr = muonCopy.fCentr;
164 0 : fPgen = muonCopy.fPgen;
165 0 : fTrackPythiaLine = muonCopy.fTrackPythiaLine;
166 0 : fTrackPDGCode = muonCopy.fTrackPDGCode;
167 0 : fOscillation = muonCopy.fOscillation;
168 0 : fNParents = muonCopy.fNParents;
169 0 : fWeight = muonCopy.fWeight;
170 :
171 0 : for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
172 0 : for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
173 0 : fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
174 0 : fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
175 : }
176 0 : for (Int_t i = 0; i < 4; i++){
177 0 : fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
178 0 : fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
179 : }
180 :
181 0 : return *this;
182 0 : }
183 :
184 : //============================================
185 :
186 : void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
187 : /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
188 0 : AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
189 0 : if (!trPar) {
190 0 : AliError("The track must contain the parameters at vertex");
191 0 : return;
192 : }
193 0 : this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
194 0 : this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
195 0 : this->SetTriggered(trackReco->GetMatchTrigger());
196 :
197 0 : Double_t xyz[3] = { trPar->GetNonBendingCoor(),
198 0 : trPar->GetBendingCoor(),
199 0 : trPar->GetZ()};
200 0 : if (zvert!=-9999) xyz[2] = zvert;
201 0 : this->SetVertex(xyz);
202 0 : }
203 :
204 : //============================================
205 : void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
206 : /// computes prec and charge from ESD track
207 0 : Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
208 0 : Double_t thetaX = muonTrack->GetThetaX();
209 0 : Double_t thetaY = muonTrack->GetThetaY();
210 0 : Double_t tanthx = TMath::Tan(thetaX);
211 0 : Double_t tanthy = TMath::Tan(thetaY);
212 0 : Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
213 0 : Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
214 0 : Double_t px = pz * tanthx;
215 0 : Double_t py = pz * tanthy;
216 0 : fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
217 0 : Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
218 0 : fPrec.SetPxPyPzE(px,py,pz,energy);
219 : // get the position
220 0 : fXYZ[0] = muonTrack->GetNonBendingCoor();
221 0 : fXYZ[1] = muonTrack->GetBendingCoor();
222 0 : if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
223 0 : else fXYZ[2] = zvert;
224 : // get the chi2 per d.o.f.
225 0 : fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
226 0 : fIsTriggered = muonTrack->GetMatchTrigger();
227 0 : }
228 :
229 : //============================================
230 : void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
231 : /// set the reconstructed 4-momentum, assuming the particle is a muon
232 0 : Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
233 0 : Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
234 0 : fPrec.SetPxPyPzE(px,py,pz,energy);
235 0 : }
236 :
237 : //============================================
238 : void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
239 : /// scans the muon history to determine parents pdg code and pythia line
240 : Int_t countP = -1;
241 0 : Int_t parents[10], parLine[10];
242 0 : Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
243 :
244 : TParticle *mother;
245 : Int_t status=-1, pdg=-1;
246 0 : while(lineM >= 0){
247 :
248 0 : mother = stack->Particle(lineM); //direct mother of rec. track
249 0 : pdg = mother->GetPdgCode();//store PDG code of first mother
250 : // break if a string, gluon, quark or diquark is found
251 0 : if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
252 0 : parents[++countP] = pdg;
253 0 : parLine[countP] = lineM;
254 0 : status = mother->GetStatusCode();//get its status code to check if oscillation occured
255 0 : if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
256 0 : lineM = mother->GetFirstMother();
257 : }
258 : //store all the fragmented parents in an array:
259 0 : for(int i = 0; i <= countP; i++){
260 0 : this->SetParentPDGCode(i,parents[countP-i]);
261 0 : this->SetParentPythiaLine(i,parLine[countP-i]);
262 : }
263 0 : fNParents = countP+1;
264 : countP = -1;
265 :
266 : //and store the lines of the string and further quarks in another array:
267 0 : while(lineM >= 0){
268 0 : mother = stack->Particle(lineM);
269 0 : pdg = mother->GetPdgCode();
270 : //now, get information before the fragmentation
271 0 : this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
272 0 : this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
273 0 : lineM = mother->GetFirstMother();
274 : }
275 :
276 : //check if in case of HF production, the string points to the correct end
277 : //and correct it in case of need:
278 : countP = 1;
279 0 : for(int par = 0; par < 4; par++){
280 0 : if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
281 : countP = par; //get the quark just before hadronisation
282 0 : break;
283 : }
284 : }
285 0 : if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
286 0 : if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
287 :
288 0 : AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
289 : this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
290 : );
291 :
292 0 : pdg = this->GetQuarkPDGCode(countP);
293 0 : Int_t line = this->GetQuarkPythiaLine(countP);
294 0 : this->ResetQuarkInfo();
295 0 : while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
296 : //must coincide with the flavour of the last fragmented mother
297 :
298 0 : pdg = stack->Particle(++line)->GetPdgCode();
299 : }
300 : //now, we have the correct string end and the correct line
301 : //continue to fill again all parents and corresponding lines
302 0 : while(line >= 0){
303 0 : mother = stack->Particle(line);//get again the mother
304 0 : pdg = mother->GetPdgCode();
305 0 : this->SetQuarkPythiaLine(countP, line);
306 0 : this->SetQuarkPDGCode(countP++, pdg);
307 0 : line = mother->GetFirstMother();
308 : }
309 0 : this->PrintInfo("h");
310 0 : }//mismatch
311 : }
312 0 : }
313 :
314 : //====================================
315 : void AliMUONTrackLight::ResetQuarkInfo(){
316 : /// resets parton information
317 0 : for(int pos = 1; pos < 4; pos++){//[0] is the string
318 0 : this->SetQuarkPDGCode(pos,-1);
319 0 : this->SetQuarkPythiaLine(pos,-1);
320 : }
321 0 : }
322 :
323 : //====================================
324 : Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
325 : /// checks if the particle is a B0
326 0 : Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
327 : Bool_t answer = kFALSE;
328 0 : for(int i = 0; i < 2; i++){
329 0 : if(TMath::Abs(intTest) == bMes0[i]){
330 : answer = kTRUE;
331 0 : break;
332 : }
333 : }
334 0 : return answer;
335 0 : }
336 : //====================================
337 : Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
338 : /// checks if mother is a resonance
339 0 : Int_t intTest = GetParentPDGCode(index);
340 : // the resonance pdg code is built this way
341 : // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
342 0 : Int_t id=intTest%100000;
343 0 : return (!((id-id%10)%110));
344 : }
345 : //====================================
346 : Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
347 : /// returns the flavour of parent idParent (idParent=0 is the oldest
348 : /// hadronized parent)
349 0 : Int_t pdg = GetParentPDGCode(idParent);
350 0 : Int_t quark = TMath::Abs(pdg/100);
351 0 : if(quark > 9) quark = quark/10;
352 0 : return quark;
353 : }
354 :
355 : //====================================
356 : void AliMUONTrackLight::PrintInfo(const Option_t* opt){
357 : /// prints information about the track:
358 : /// - "H" muon's decay history
359 : /// - "K" muon kinematics
360 : /// - "A" all variables
361 0 : TString options(opt);
362 0 : options.ToUpper();
363 :
364 0 : if(options.Contains("H") || options.Contains("A")){ //muon decay history
365 0 : char *name= new char[100];
366 0 : TString pdg = "", line = "";
367 0 : for(int i = 3; i >= 0; i--){
368 0 : if(this->GetQuarkPythiaLine(i)>= 0){
369 0 : snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
370 0 : line += name;
371 0 : snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
372 0 : pdg += name;
373 : }
374 : }
375 0 : for(int i = 0; i < fNParents; i++){
376 0 : if(this->GetParentPythiaLine(i)>= 0){
377 0 : snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
378 0 : line += name;
379 0 : snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
380 0 : pdg += name;
381 : }
382 : }
383 0 : snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
384 0 : snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
385 :
386 0 : printf("\nmuon's decay history:\n");
387 0 : printf(" PDG: %s\n", pdg.Data());
388 0 : printf("line: %s\n", line.Data());
389 0 : }
390 0 : if(options.Contains("K") || options.Contains("A")){ //muon kinematic
391 :
392 0 : Int_t charge = this->GetCharge();
393 0 : Double_t *vtx = this->GetVertex();
394 0 : TLorentzVector momRec = this->GetPRec();
395 0 : TLorentzVector momGen = this->GetPGen();
396 0 : printf("the track's charge is %d\n", charge);
397 0 : printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
398 0 : printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
399 0 : printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
400 0 : printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
401 0 : momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
402 0 : momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
403 0 : }
404 0 : }
405 : //====================================
406 : Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
407 : /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
408 0 : Int_t pdg = this->GetParentPDGCode(idparent);
409 0 : if (TMath::Abs(pdg)==211 || //pi+
410 0 : TMath::Abs(pdg)==321 || //K+
411 0 : TMath::Abs(pdg)==213 || //rho+
412 0 : TMath::Abs(pdg)==311 || //K0
413 0 : TMath::Abs(pdg)==313 || //K*0
414 0 : TMath::Abs(pdg)==323 //K*+
415 : ) {
416 0 : return kTRUE;
417 : }
418 0 : else return kFALSE;
419 0 : }
420 : //====================================
421 : Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
422 : /// check if the provided pdg code corresponds to a diquark
423 0 : pdg = TMath::Abs(pdg);
424 0 : if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
425 0 : else return kFALSE;
426 0 : }
427 : //====================================
428 : void AliMUONTrackLight::SetParentPDGCode(Int_t index, Int_t pdg) {
429 : /// Set hadronised parents and grandparents
430 0 : if ( index < fgkNParentsMax ) {
431 0 : fParentPDGCode[index] = pdg;
432 0 : } else {
433 0 : AliErrorStream() << "Index outside the array size." << std::endl;
434 : }
435 0 : }
436 : //====================================
437 : void AliMUONTrackLight::SetParentPythiaLine(Int_t index, Int_t line) {
438 : /// Set line of Pythia output for hadronised parents & grandparents
439 0 : if ( index < fgkNParentsMax ) {
440 0 : fParentPythiaLine[index] = line;
441 0 : } else {
442 0 : AliErrorStream() << "Index outside the array size." << std::endl;
443 : }
444 0 : }
445 : //====================================
446 : void AliMUONTrackLight::SetQuarkPDGCode(Int_t index, Int_t pdg){
447 : /// Set pdg of the string [0], quarks/gluons [1,2], sometimes proton [3]
448 0 : if ( index < 4 ) {
449 0 : fQuarkPDGCode[index] = pdg;
450 0 : } else {
451 0 : AliErrorStream() << "Index outside the array size." << std::endl;
452 : }
453 0 : }
454 : //====================================
455 : void AliMUONTrackLight::SetQuarkPythiaLine(Int_t index, Int_t line){
456 : /// Set line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
457 0 : if ( index < 4 ) {
458 0 : fQuarkPythiaLine[index] = line;
459 0 : } else {
460 0 : AliErrorStream() << "Index outside the array size." << std::endl;
461 : }
462 0 : }
|
