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 #include "AliMUONTrackExtrap.h"
  42
  43 #include "TDatabasePDG.h"
  44 #include "TClonesArray.h"
  45 #include "TParticle.h"
  46 #include "TString.h"
  47
  48 ClassImp(AliMUONTrackLight)
  49
  50 //===================================================================
  51
  52 AliMUONTrackLight::AliMUONTrackLight()
  53   : TObject(),
  54     fPrec(),
  55     fIsTriggered(kFALSE),
  56     fCharge(-999),
  57     fChi2(-1),
  58     fCentr(-1),
  59     fPgen(),
  60     fTrackPythiaLine(-999),
  61     fTrackPDGCode(-999),
  62     fOscillation(kFALSE),
  63     fNParents(0),
  64     fWeight(1)
  65 {
  66   /// default constructor
  67   fPgen.SetPxPyPzE(0.,0.,0.,0.);
  68   fPrec.SetPxPyPzE(0.,0.,0.,0.);
  69   for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
  70   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  71     fParentPDGCode[npar] = -1;
  72     fParentPythiaLine[npar] = -1;
  73   }
  74   for (Int_t i = 0; i < 4; i++){
  75     fQuarkPDGCode[i] = -1;
  76     fQuarkPythiaLine[i] = -1;
  77   }
  78 }
  79
  80 //============================================
  81 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
  82   : TObject(muonCopy),
  83     fPrec(muonCopy.fPrec),
  84     fIsTriggered(muonCopy.fIsTriggered),
  85     fCharge(muonCopy.fCharge),
  86     fChi2(muonCopy.fChi2),
  87     fCentr(muonCopy.fCentr),
  88     fPgen(muonCopy.fPgen),
  89     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
  90     fTrackPDGCode(muonCopy.fTrackPDGCode),
  91     fOscillation(muonCopy.fOscillation),
  92     fNParents(muonCopy.fNParents),
  93     fWeight(muonCopy.fWeight)
  94 {
  95   /// copy constructor
  96   for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
  97   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  98     fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
  99     fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
 100   }
 101   for (Int_t i = 0; i < 4; i++){
 102     fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
 103     fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
 104   }
 105 }
 106
 107 //============================================
 108 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
 109   : TObject(),
 110     fPrec(),
 111     fIsTriggered(kFALSE),
 112     fCharge(-999),
 113     fChi2(-1),
 114     fCentr(-1),
 115     fPgen(),
 116     fTrackPythiaLine(-999),
 117     fTrackPDGCode(-999),
 118     fOscillation(kFALSE),
 119     fNParents(0),
 120     fWeight(1)
 121 {
 122   /// constructor
 123   //AliMUONTrackLight();
 124   FillFromESD(muonTrack);
 125 }
 126
 127 //============================================
 128 AliMUONTrackLight::~AliMUONTrackLight()
 129 {
 130 /// Destructor
 131 }
 132
 133 //============================================
 134
 135 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
 136   /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
 137   AliMUONTrackParam trPar(*((AliMUONTrackParam*) (trackReco->GetTrackParamAtHit()->First())));
 138   //  AliMUONTrackParam *trPar  = trackReco->GetTrackParamAtVertex();
 139   AliMUONTrackExtrap::ExtrapToVertex(&trPar,0.,0.,0.);
 140   this->SetCharge(Int_t(TMath::Sign(1.,trPar.GetInverseBendingMomentum())));
 141   this->SetPxPyPz(trPar.Px(),trPar.Py(), trPar.Pz());
 142   this->SetTriggered(trackReco->GetMatchTrigger());
 143
 144   Double_t xyz[3] = { trPar.GetNonBendingCoor(),
 145                       trPar.GetBendingCoor(),
 146                       zvert};
 147   this->SetVertex(xyz);
 148 }
 149
 150 //============================================
 151 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
 152   /// computes prec and charge from ESD track
 153   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 154   Double_t thetaX = muonTrack->GetThetaX();
 155   Double_t thetaY = muonTrack->GetThetaY();
 156   Double_t tanthx = TMath::Tan(thetaX);
 157   Double_t tanthy = TMath::Tan(thetaY);
 158   Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 159   Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
 160   Double_t px     = pz * tanthx;
 161   Double_t py     = pz * tanthy;
 162   fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 163   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 164   fPrec.SetPxPyPzE(px,py,pz,energy);
 165   // get the position
 166   fXYZ[0] = muonTrack->GetNonBendingCoor();
 167   fXYZ[1] = muonTrack->GetBendingCoor();
 168   if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
 169   else fXYZ[2] = zvert;
 170   // get the chi2 per d.o.f.
 171   fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
 172   fIsTriggered = muonTrack->GetMatchTrigger();
 173 }
 174
 175 //============================================
 176 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
 177   /// set the reconstructed 4-momentum, assuming the particle is a muon
 178   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 179   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 180   fPrec.SetPxPyPzE(px,py,pz,energy);
 181 }
 182
 183 //============================================
 184 TParticle* AliMUONTrackLight::FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader){
 185   /// find the MC particle that corresponds to a given rec track
 186   TParticle *part = 0;
 187   const Double_t kSigma2Cut = 16;  // 4 sigmas cut, kSigma2Cut = 4*4
 188   Int_t nTrackRef = trackRefArray->GetEntriesFast();
 189   Int_t compPart = 0;
 190   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 191     AliMUONTrack *trackRef = (AliMUONTrack *)trackRefArray->At(iref);
 192     // check if trackRef is compatible with trackReco:
 193     //routine returns for each chamber a yes/no information if the
 194     //hit of rec. track and hit of referenced track are compatible
 195     Bool_t *compTrack = trackRef->CompatibleTrack(trackReco,kSigma2Cut);
 196     Int_t iTrack = this->TrackCheck(compTrack); //returns number of validated conditions
 197     if (iTrack==4) {
 198       compPart++;
 199       Int_t trackID = trackRef->GetTrackID();
 200       this->SetTrackPythiaLine(trackID);
 201       part = ((AliStack *)(((AliHeader *) runLoader->GetHeader())->Stack()))->Particle(trackID);
 202       fTrackPDGCode = part->GetPdgCode();
 203     }
 204   }
 205   if (compPart>1) {
 206     printf ("<AliMUONTrackLight::FindRefTrack> ERROR: more than one particle compatible to the reconstructed track.\n");
 207     Int_t i=0, j=1/i;
 208     printf ("j=%d \n",j);
 209   }
 210   return part;
 211 }
 212
 213 //============================================
 214 Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
 215   /// Apply reconstruction requirements
 216   /// Return number of validated conditions
 217   /// If all the tests are verified then TrackCheck = 4 (good track)
 218   Int_t iTrack = 0;
 219   Int_t hitsInLastStations = 0;
 220
 221   // apply reconstruction requirements
 222   if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
 223   if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
 224   if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
 225   for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
 226     if (compTrack[ch]) hitsInLastStations++;
 227   }
 228   if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
 229   return iTrack;
 230 }
 231
 232 //============================================
 233 void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part){
 234   /// scans the muon history to determine parents pdg code and pythia line
 235   // kept for backward compatibility
 236   // see the overloaded method FillMuonHistory(AliStack *stack, TParticle *part)
 237   AliStack *stack = runLoader->GetHeader()->Stack();
 238   FillMuonHistory(stack,part);
 239 }
 240
 241 //============================================
 242 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
 243   /// scans the muon history to determine parents pdg code and pythia line
 244   Int_t countP = -1;
 245   Int_t parents[10], parLine[10];
 246   Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
 247
 248   TParticle *mother;
 249   Int_t status=-1, pdg=-1;
 250   while(lineM >= 0){
 251
 252     mother = stack->Particle(lineM); //direct mother of rec. track
 253     pdg = mother->GetPdgCode();//store PDG code of first mother
 254     // break if a string, gluon, quark or diquark is found
 255     if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
 256     parents[++countP] = pdg;
 257     parLine[countP] = lineM;
 258     status = mother->GetStatusCode();//get its status code to check if oscillation occured
 259     if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
 260     lineM = mother->GetFirstMother();
 261   }
 262   //store all the fragmented parents in an array:
 263   for(int i = 0; i <= countP; i++){
 264     this->SetParentPDGCode(i,parents[countP-i]);
 265     this->SetParentPythiaLine(i,parLine[countP-i]);
 266   }
 267   fNParents = countP+1;
 268   countP = -1;
 269
 270   //and store the lines of the string and further quarks in another array:
 271   while(lineM >= 0){
 272     mother = stack->Particle(lineM);
 273     pdg = mother->GetPdgCode();
 274     //now, get information before the fragmentation
 275     this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
 276     this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
 277     lineM = mother->GetFirstMother();
 278   }
 279
 280   //check if in case of HF production, the string points to the correct end
 281   //and correct it in case of need:
 282   countP = 1;
 283   for(int par = 0; par < 4; par++){
 284     if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
 285       countP = par; //get the quark just before hadronisation
 286       break;
 287     }
 288   }
 289   if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
 290     if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
 291
 292       printf("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n", this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)));
 293       Int_t pdg = this->GetQuarkPDGCode(countP), line = this->GetQuarkPythiaLine(countP);
 294       this->ResetQuarkInfo();
 295       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         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       while(line >= 0){
 303         mother = stack->Particle(line);//get again the mother
 304         pdg = mother->GetPdgCode();
 305         this->SetQuarkPythiaLine(countP, line);
 306         this->SetQuarkPDGCode(countP++, pdg);
 307         line = mother->GetFirstMother();
 308       }
 309       this->PrintInfo("h");
 310     }//mismatch
 311   }
 312 }
 313
 314 //====================================
 315 void AliMUONTrackLight::ResetQuarkInfo(){
 316   /// resets parton information
 317   for(int pos = 1; pos < 4; pos++){//[0] is the string
 318     this->SetQuarkPDGCode(pos,-1);
 319     this->SetQuarkPythiaLine(pos,-1);
 320   }
 321 }
 322
 323 //====================================
 324 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
 325   /// checks if the particle is a B0
 326   Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
 327   Bool_t answer = kFALSE;
 328   for(int i = 0; i < 2; i++){
 329     if(TMath::Abs(intTest) == bMes0[i]){
 330       answer = kTRUE;
 331       break;
 332     }
 333   }
 334   return answer;
 335 }
 336 //====================================
 337 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
 338   /// checks if mother is a resonance
 339   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   Int_t id=intTest%100000;
 343   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   Int_t pdg = GetParentPDGCode(idParent);
 350   Int_t quark = TMath::Abs(pdg/100);
 351   if(quark > 9) quark = quark/10;
 352   return quark;
 353 }
 354
 355 //====================================
 356 void AliMUONTrackLight::PrintInfo(Option_t* opt){
 357   /// prints information about the track:
 358   /// - "H" muon's decay history
 359   /// - "K" muon kinematics
 360   /// - "A" all variables
 361   TString options(opt);
 362   options.ToUpper();
 363
 364   if(options.Contains("H") || options.Contains("A")){ //muon decay history
 365     char *name= new char[100];
 366     TString pdg = "", line = "";
 367     for(int i = 3; i >= 0; i--){
 368       if(this->GetQuarkPythiaLine(i)>= 0){
 369         sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
 370         line += name;
 371         sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
 372         pdg += name;
 373       }
 374     }
 375     for(int i = 0; i < fNParents; i++){
 376       if(this->GetParentPythiaLine(i)>= 0){
 377         sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
 378         line += name;
 379         sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
 380         pdg += name;
 381       }
 382     }
 383     sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
 384     sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
 385
 386     printf("\nmuon's decay history:\n");
 387     printf(" PDG: %s\n", pdg.Data());
 388     printf("line: %s\n", line.Data());
 389   }
 390   if(options.Contains("K") || options.Contains("A")){ //muon kinematic
 391
 392     Int_t charge = this->GetCharge();
 393     Double_t *vtx = this->GetVertex();
 394     TLorentzVector momRec = this->GetPRec();
 395     TLorentzVector momGen = this->GetPGen();
 396     printf("the track's charge is %d\n", charge);
 397     printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 398     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 399     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 400     printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
 401            momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 402            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 403   }
 404 }
 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   Int_t pdg = this->GetParentPDGCode(idparent);
 409   if (TMath::Abs(pdg)==211 || //pi+
 410       TMath::Abs(pdg)==321 || //K+
 411       TMath::Abs(pdg)==213 || //rho+
 412       TMath::Abs(pdg)==311 || //K0
 413       TMath::Abs(pdg)==313 || //K*0
 414       TMath::Abs(pdg)==323    //K*+
 415       ) {
 416     return kTRUE;
 417   }
 418   else return kFALSE;
 419 }
 420 //====================================
 421 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg){
 422   /// check if the provided pdg code corresponds to a diquark
 423   pdg = TMath::Abs(pdg);
 424   if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
 425   else return kFALSE;
 426 }