MUON/AliMUONTrackLight.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 // 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 ClassImp(AliMUONTrackLight)
  54
  55 //===================================================================
  56
  57 AliMUONTrackLight::AliMUONTrackLight()
  58   : TObject(),
  59     fPrec(),
  60     fIsTriggered(kFALSE),
  61     fCharge(-999),
  62     fChi2(-1),
  63     fCentr(-1),
  64     fPgen(),
  65     fTrackPythiaLine(-999),
  66     fTrackPDGCode(-999),
  67     fOscillation(kFALSE),
  68     fNParents(0),
  69     fWeight(1)
  70 {
  71   /// default constructor
  72   fPgen.SetPxPyPzE(0.,0.,0.,0.);
  73   fPrec.SetPxPyPzE(0.,0.,0.,0.);
  74   for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
  75   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  76     fParentPDGCode[npar] = -1;
  77     fParentPythiaLine[npar] = -1;
  78   }
  79   for (Int_t i = 0; i < 4; i++){
  80     fQuarkPDGCode[i] = -1;
  81     fQuarkPythiaLine[i] = -1;
  82   }
  83 }
  84
  85 //============================================
  86 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
  87   : TObject(muonCopy),
  88     fPrec(muonCopy.fPrec),
  89     fIsTriggered(muonCopy.fIsTriggered),
  90     fCharge(muonCopy.fCharge),
  91     fChi2(muonCopy.fChi2),
  92     fCentr(muonCopy.fCentr),
  93     fPgen(muonCopy.fPgen),
  94     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
  95     fTrackPDGCode(muonCopy.fTrackPDGCode),
  96     fOscillation(muonCopy.fOscillation),
  97     fNParents(muonCopy.fNParents),
  98     fWeight(muonCopy.fWeight)
  99 {
 100   /// copy constructor
 101   for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
 102   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
 103     fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
 104     fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
 105   }
 106   for (Int_t i = 0; i < 4; i++){
 107     fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
 108     fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
 109   }
 110 }
 111
 112 //============================================
 113 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
 114   : TObject(),
 115     fPrec(),
 116     fIsTriggered(kFALSE),
 117     fCharge(-999),
 118     fChi2(-1),
 119     fCentr(-1),
 120     fPgen(),
 121     fTrackPythiaLine(-999),
 122     fTrackPDGCode(-999),
 123     fOscillation(kFALSE),
 124     fNParents(0),
 125     fWeight(1)
 126 {
 127   /// constructor
 128   fPgen.SetPxPyPzE(0.,0.,0.,0.);
 129   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
 130     fParentPDGCode[npar] = -1;
 131     fParentPythiaLine[npar] = -1;
 132   }
 133   for (Int_t i = 0; i < 4; i++){
 134     fQuarkPDGCode[i] = -1;
 135     fQuarkPythiaLine[i] = -1;
 136   }
 137   FillFromESD(muonTrack);
 138 }
 139
 140 //============================================
 141 AliMUONTrackLight::~AliMUONTrackLight()
 142 {
 143 /// Destructor
 144 }
 145
 146 //============================================
 147
 148 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
 149   /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
 150   AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
 151   if (!trPar) {
 152     AliError("The track must contain the parameters at vertex");
 153     return;
 154   }
 155   this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
 156   this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
 157   this->SetTriggered(trackReco->GetMatchTrigger());
 158
 159   Double_t xyz[3] = { trPar->GetNonBendingCoor(),
 160                       trPar->GetBendingCoor(),
 161                       trPar->GetZ()};
 162   if (zvert!=-9999) xyz[2] = zvert;
 163   this->SetVertex(xyz);
 164 }
 165
 166 //============================================
 167 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
 168   /// computes prec and charge from ESD track
 169   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 170   Double_t thetaX = muonTrack->GetThetaX();
 171   Double_t thetaY = muonTrack->GetThetaY();
 172   Double_t tanthx = TMath::Tan(thetaX);
 173   Double_t tanthy = TMath::Tan(thetaY);
 174   Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 175   Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
 176   Double_t px     = pz * tanthx;
 177   Double_t py     = pz * tanthy;
 178   fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 179   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 180   fPrec.SetPxPyPzE(px,py,pz,energy);
 181   // get the position
 182   fXYZ[0] = muonTrack->GetNonBendingCoor();
 183   fXYZ[1] = muonTrack->GetBendingCoor();
 184   if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
 185   else fXYZ[2] = zvert;
 186   // get the chi2 per d.o.f.
 187   fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
 188   fIsTriggered = muonTrack->GetMatchTrigger();
 189 }
 190
 191 //============================================
 192 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
 193   /// set the reconstructed 4-momentum, assuming the particle is a muon
 194   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 195   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 196   fPrec.SetPxPyPzE(px,py,pz,energy);
 197 }
 198
 199 //============================================
 200 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
 201   /// scans the muon history to determine parents pdg code and pythia line
 202   Int_t countP = -1;
 203   Int_t parents[10], parLine[10];
 204   Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
 205
 206   TParticle *mother;
 207   Int_t status=-1, pdg=-1;
 208   while(lineM >= 0){
 209
 210     mother = stack->Particle(lineM); //direct mother of rec. track
 211     pdg = mother->GetPdgCode();//store PDG code of first mother
 212     // break if a string, gluon, quark or diquark is found
 213     if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
 214     parents[++countP] = pdg;
 215     parLine[countP] = lineM;
 216     status = mother->GetStatusCode();//get its status code to check if oscillation occured
 217     if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
 218     lineM = mother->GetFirstMother();
 219   }
 220   //store all the fragmented parents in an array:
 221   for(int i = 0; i <= countP; i++){
 222     this->SetParentPDGCode(i,parents[countP-i]);
 223     this->SetParentPythiaLine(i,parLine[countP-i]);
 224   }
 225   fNParents = countP+1;
 226   countP = -1;
 227
 228   //and store the lines of the string and further quarks in another array:
 229   while(lineM >= 0){
 230     mother = stack->Particle(lineM);
 231     pdg = mother->GetPdgCode();
 232     //now, get information before the fragmentation
 233     this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
 234     this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
 235     lineM = mother->GetFirstMother();
 236   }
 237
 238   //check if in case of HF production, the string points to the correct end
 239   //and correct it in case of need:
 240   countP = 1;
 241   for(int par = 0; par < 4; par++){
 242     if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
 243       countP = par; //get the quark just before hadronisation
 244       break;
 245     }
 246   }
 247   if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
 248     if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
 249
 250       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 251           this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
 252         );
 253
 254       pdg = this->GetQuarkPDGCode(countP);
 255       Int_t line = this->GetQuarkPythiaLine(countP);
 256       this->ResetQuarkInfo();
 257       while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
 258                                                         //must coincide with the flavour of the last fragmented mother
 259
 260         pdg = stack->Particle(++line)->GetPdgCode();
 261       }
 262       //now, we have the correct string end and the correct line
 263       //continue to fill again all parents and corresponding lines
 264       while(line >= 0){
 265         mother = stack->Particle(line);//get again the mother
 266         pdg = mother->GetPdgCode();
 267         this->SetQuarkPythiaLine(countP, line);
 268         this->SetQuarkPDGCode(countP++, pdg);
 269         line = mother->GetFirstMother();
 270       }
 271       this->PrintInfo("h");
 272     }//mismatch
 273   }
 274 }
 275
 276 //====================================
 277 void AliMUONTrackLight::ResetQuarkInfo(){
 278   /// resets parton information
 279   for(int pos = 1; pos < 4; pos++){//[0] is the string
 280     this->SetQuarkPDGCode(pos,-1);
 281     this->SetQuarkPythiaLine(pos,-1);
 282   }
 283 }
 284
 285 //====================================
 286 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
 287   /// checks if the particle is a B0
 288   Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
 289   Bool_t answer = kFALSE;
 290   for(int i = 0; i < 2; i++){
 291     if(TMath::Abs(intTest) == bMes0[i]){
 292       answer = kTRUE;
 293       break;
 294     }
 295   }
 296   return answer;
 297 }
 298 //====================================
 299 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
 300   /// checks if mother is a resonance
 301   Int_t intTest = GetParentPDGCode(index);
 302   // the resonance pdg code is built this way
 303   // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
 304   Int_t id=intTest%100000;
 305   return (!((id-id%10)%110));
 306 }
 307 //====================================
 308 Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
 309   /// returns the flavour of parent idParent (idParent=0 is the oldest
 310   /// hadronized parent)
 311   Int_t pdg = GetParentPDGCode(idParent);
 312   Int_t quark = TMath::Abs(pdg/100);
 313   if(quark > 9) quark = quark/10;
 314   return quark;
 315 }
 316
 317 //====================================
 318 void AliMUONTrackLight::PrintInfo(Option_t* opt){
 319   /// prints information about the track:
 320   /// - "H" muon's decay history
 321   /// - "K" muon kinematics
 322   /// - "A" all variables
 323   TString options(opt);
 324   options.ToUpper();
 325
 326   if(options.Contains("H") || options.Contains("A")){ //muon decay history
 327     char *name= new char[100];
 328     TString pdg = "", line = "";
 329     for(int i = 3; i >= 0; i--){
 330       if(this->GetQuarkPythiaLine(i)>= 0){
 331         sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
 332         line += name;
 333         sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
 334         pdg += name;
 335       }
 336     }
 337     for(int i = 0; i < fNParents; i++){
 338       if(this->GetParentPythiaLine(i)>= 0){
 339         sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
 340         line += name;
 341         sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
 342         pdg += name;
 343       }
 344     }
 345     sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
 346     sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
 347
 348     printf("\nmuon's decay history:\n");
 349     printf(" PDG: %s\n", pdg.Data());
 350     printf("line: %s\n", line.Data());
 351   }
 352   if(options.Contains("K") || options.Contains("A")){ //muon kinematic
 353
 354     Int_t charge = this->GetCharge();
 355     Double_t *vtx = this->GetVertex();
 356     TLorentzVector momRec = this->GetPRec();
 357     TLorentzVector momGen = this->GetPGen();
 358     printf("the track's charge is %d\n", charge);
 359     printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 360     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 361     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 362     printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
 363            momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 364            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 365   }
 366 }
 367 //====================================
 368 Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
 369   /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
 370   Int_t pdg = this->GetParentPDGCode(idparent);
 371   if (TMath::Abs(pdg)==211 || //pi+
 372       TMath::Abs(pdg)==321 || //K+
 373       TMath::Abs(pdg)==213 || //rho+
 374       TMath::Abs(pdg)==311 || //K0
 375       TMath::Abs(pdg)==313 || //K*0
 376       TMath::Abs(pdg)==323    //K*+
 377       ) {
 378     return kTRUE;
 379   }
 380   else return kFALSE;
 381 }
 382 //====================================
 383 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg){
 384   /// check if the provided pdg code corresponds to a diquark
 385   pdg = TMath::Abs(pdg);
 386   if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
 387   else return kFALSE;
 388 }