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 //This class was prepared by INFN Cagliari, July 2006
  20 //(authors: H.Woehri, A.de Falco)
  21 //
  22 //
  23 // Compact information for the muon generated tracks in the MUON arm
  24 // useful at the last stage of the analysis chain
  25 // provides a link between the reconstructed track and the generated particle
  26 // stores kinematical information at gen. and rec. level and
  27 // the decay history of the muon, allowing the identification of the
  28 // mother process
  29 //
  30 // To be used together with AliMUONPairLight
  31 //===================================================================
  32
  33 #include "AliMUONTrackLight.h"
  34 #include "AliMUONTrack.h"
  35 #include "AliMUONConstants.h"
  36
  37 #include "AliESDMuonTrack.h"
  38 #include "AliRunLoader.h"
  39 #include "AliStack.h"
  40 #include "AliHeader.h"
  41
  42 #include "TDatabasePDG.h"
  43 #include "TClonesArray.h"
  44 #include "TParticle.h"
  45 #include "TString.h"
  46
  47 ClassImp(AliMUONTrackLight)
  48
  49 //===================================================================
  50
  51 AliMUONTrackLight::AliMUONTrackLight()
  52   : TObject(),
  53     fPrec(),
  54     fIsTriggered(kFALSE),
  55     fCharge(-999),
  56     fCentr(-1),
  57     fPgen(),
  58     fTrackPythiaLine(-999),
  59     fTrackPDGCode(-999),
  60     fOscillation(kFALSE),
  61     fNParents(0)
  62 {
  63   /// default constructor
  64   fPgen.SetPxPyPzE(0.,0.,0.,0.);
  65   fPrec.SetPxPyPzE(0.,0.,0.,0.);
  66   for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
  67   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  68     fParentPDGCode[npar] = -1;
  69     fParentPythiaLine[npar] = -1;
  70   }
  71   for (Int_t i = 0; i < 4; i++){
  72     fQuarkPDGCode[i] = -1;
  73     fQuarkPythiaLine[i] = -1;
  74   }
  75 }
  76
  77 //============================================
  78 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
  79   : TObject(muonCopy),
  80     fPrec(muonCopy.fPrec),
  81     fIsTriggered(muonCopy.fIsTriggered),
  82     fCharge(muonCopy.fCharge),
  83     fCentr(muonCopy.fCentr),
  84     fPgen(muonCopy.fPgen),
  85     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
  86     fTrackPDGCode(muonCopy.fTrackPDGCode),
  87     fOscillation(muonCopy.fOscillation),
  88     fNParents(muonCopy.fNParents)
  89 {
  90   /// copy constructor
  91   for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
  92   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  93     fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
  94     fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
  95   }
  96   for (Int_t i = 0; i < 4; i++){
  97     fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
  98     fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
  99   }
 100 }
 101
 102 //============================================
 103 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
 104   : TObject(),
 105     fPrec(),
 106     fIsTriggered(kFALSE),
 107     fCharge(-999),
 108     fCentr(-1),
 109     fPgen(),
 110     fTrackPythiaLine(-999),
 111     fTrackPDGCode(-999),
 112     fOscillation(kFALSE),
 113     fNParents(0)
 114 {
 115   /// constructor
 116   //AliMUONTrackLight();
 117   ComputePRecAndChargeFromESD(muonTrack);
 118 }
 119
 120 //============================================
 121 void AliMUONTrackLight::ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack){
 122   /// computes prec and charge from ESD track
 123   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 124   Double_t thetaX = muonTrack->GetThetaX();
 125   Double_t thetaY = muonTrack->GetThetaY();
 126   Double_t tanthx = TMath::Tan(thetaX);
 127   Double_t tanthy = TMath::Tan(thetaY);
 128   Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 129   Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
 130   Double_t px     = pz * tanthx;
 131   Double_t py     = pz * tanthy;
 132   fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 133   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 134   fPrec.SetPxPyPzE(px,py,pz,energy);
 135 }
 136
 137 //============================================
 138 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
 139   /// set the reconstructed 4-momentum, assuming the particle is a muon
 140   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 141   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 142   fPrec.SetPxPyPzE(px,py,pz,energy);
 143 }
 144
 145 //============================================
 146 TParticle* AliMUONTrackLight::FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader){
 147   /// find the MC particle that corresponds to a given rec track
 148   TParticle *part = 0;
 149   const Double_t kSigma2Cut = 16;  // 4 sigmas cut, kSigma2Cut = 4*4
 150   Int_t nTrackRef = trackRefArray->GetEntriesFast();
 151   Int_t compPart = 0;
 152   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 153     AliMUONTrack *trackRef = (AliMUONTrack *)trackRefArray->At(iref);
 154     // check if trackRef is compatible with trackReco:
 155     //routine returns for each chamber a yes/no information if the
 156     //hit of rec. track and hit of referenced track are compatible
 157     Bool_t *compTrack = trackRef->CompatibleTrack(trackReco,kSigma2Cut);
 158     Int_t iTrack = this->TrackCheck(compTrack); //returns number of validated conditions
 159     if (iTrack==4) {
 160       compPart++;
 161       Int_t trackID = trackRef->GetTrackID();
 162       this->SetTrackPythiaLine(trackID);
 163       part = ((AliStack *)(((AliHeader *) runLoader->GetHeader())->Stack()))->Particle(trackID);
 164       fTrackPDGCode = part->GetPdgCode();
 165     }
 166   }
 167   if (compPart>1) {
 168     printf ("<AliMUONTrackLight::FindRefTrack> ERROR: more than one particle compatible to the reconstructed track.\n");
 169     Int_t i=0, j=1/i;
 170     printf ("j=%d \n",j);
 171   }
 172   return part;
 173 }
 174
 175 //============================================
 176 Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
 177   /// Apply reconstruction requirements
 178   /// Return number of validated conditions
 179   /// If all the tests are verified then TrackCheck = 4 (good track)
 180   Int_t iTrack = 0;
 181   Int_t hitsInLastStations = 0;
 182
 183   // apply reconstruction requirements
 184   if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
 185   if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
 186   if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
 187   for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
 188     if (compTrack[ch]) hitsInLastStations++;
 189   }
 190   if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
 191   return iTrack;
 192 }
 193
 194 //============================================
 195 void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part){
 196   /// scans the muon history to determine parents pdg code and pythia line
 197   Int_t countP = -1;
 198   Int_t parents[10], parLine[10];
 199   Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
 200
 201   AliStack *stack = runLoader->GetHeader()->Stack();
 202   TParticle *mother;
 203   Int_t status=-1, pdg=-1;
 204   while(lineM >= 0){
 205
 206     mother = stack->Particle(lineM); //direct mother of rec. track
 207     pdg = mother->GetPdgCode();//store PDG code of first mother
 208     if(pdg == 92) break;  // break if a string is found
 209     parents[++countP] = pdg;
 210     parLine[countP] = lineM;
 211     status = mother->GetStatusCode();//get its status code to check if oscillation occured
 212     if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
 213     lineM = mother->GetFirstMother();
 214   }
 215   //store all the fragmented parents in an array:
 216   for(int i = 0; i <= countP; i++){
 217     this->SetParentPDGCode(i,parents[countP-i]);
 218     this->SetParentPythiaLine(i,parLine[countP-i]);
 219   }
 220
 221   fNParents = countP+1;
 222   countP = -1;
 223   //and store the lines of the string and further quarks in another array:
 224   while(lineM >= 0){
 225     mother = stack->Particle(lineM);
 226     pdg = mother->GetPdgCode();
 227     //now, get information before the fragmentation
 228     this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
 229     this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
 230     lineM = mother->GetFirstMother();
 231   }
 232
 233   //check if in case of HF production, the string points to the correct end
 234   //and correct it in case of need:
 235   countP = 1;
 236   if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
 237     if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
 238
 239       Int_t pdg = this->GetQuarkPDGCode(1), line = this->GetQuarkPythiaLine(1);
 240       this->ResetQuarkInfo();
 241       while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
 242                                                         //must coincide with the flavour of the last fragmented mother
 243
 244         pdg = stack->Particle(++line)->GetPdgCode();
 245       }
 246       //now, we have the correct string end and the correct line
 247       //continue to fill again all parents and corresponding lines
 248       while(line >= 0){
 249         mother = stack->Particle(line);//get again the mother
 250         pdg = mother->GetPdgCode();
 251         this->SetQuarkPythiaLine(countP, line);
 252         this->SetQuarkPDGCode(countP++, pdg);
 253         line = mother->GetFirstMother();
 254       }
 255     }
 256   }
 257 }
 258 //====================================
 259 void AliMUONTrackLight::ResetQuarkInfo(){
 260   /// resets parton information
 261   for(int pos = 1; pos < 4; pos++){//[0] is the string
 262     this->SetQuarkPDGCode(pos,-1);
 263     this->SetQuarkPythiaLine(pos,-1);
 264   }
 265 }
 266
 267 //====================================
 268 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
 269   /// checks if the particle is a B0
 270   Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
 271   Bool_t answer = kFALSE;
 272   for(int i = 0; i < 2; i++){
 273     if(TMath::Abs(intTest) == bMes0[i]){
 274       answer = kTRUE;
 275       break;
 276     }
 277   }
 278   return answer;
 279 }
 280 //====================================
 281 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
 282   /// checks if mother is a resonance
 283   Int_t intTest = GetParentPDGCode(index);
 284   // the resonance pdg code is built this way
 285   // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
 286   Int_t id=intTest%100000;
 287   return (!((id-id%10)%110));
 288 }
 289 //====================================
 290 Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
 291   /// returns the flavour of parent idParent (idParent=0 is the oldest
 292   /// hadronized parent)
 293   Int_t pdg = GetParentPDGCode(idParent);
 294   Int_t quark = TMath::Abs(pdg/100);
 295   if(quark > 9) quark = quark/10;
 296   return quark;
 297 }
 298
 299 //====================================
 300 void AliMUONTrackLight::PrintInfo(Option_t* opt){
 301   /// prints information about the track:
 302   /// - "H" muon's decay history
 303   /// - "K" muon kinematics
 304   /// - "A" all variables
 305   TString options(opt);
 306   options.ToUpper();
 307
 308   if(options.Contains("H") || options.Contains("A")){ //muon decay history
 309     char *name= new char[100];
 310     TString pdg = "", line = "";
 311     for(int i = 3; i >= 0; i--){
 312       if(this->GetQuarkPythiaLine(i)>= 0){
 313         sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
 314         line += name;
 315         sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
 316         pdg += name;
 317       }
 318     }
 319     for(int i = 0; i < fNParents; i++){
 320       if(this->GetParentPythiaLine(i)>= 0){
 321         sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
 322         line += name;
 323         sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
 324         pdg += name;
 325       }
 326     }
 327     sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
 328     sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
 329
 330     printf("\nmuon's decay history:\n");
 331     printf(" PDG: %s\n", pdg.Data());
 332     printf("line: %s\n", line.Data());
 333   }
 334   if(options.Contains("K") || options.Contains("A")){ //muon kinematic
 335
 336     Int_t charge = this->GetCharge();
 337     Double_t *vtx = this->GetVertex();
 338     TLorentzVector momRec = this->GetPRec();
 339     TLorentzVector momGen = this->GetPGen();
 340     printf("the track's charge is %d\n", charge);
 341     printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 342     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 343     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 344     printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
 345            momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 346            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 347   }
 348 }
 349
 350