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 "AliMUONTrackExtrap.h"
  44
  45 #include "AliESDMuonTrack.h"
  46 #include "AliRunLoader.h"
  47 #include "AliStack.h"
  48 #include "AliHeader.h"
  49
  50 #include "TDatabasePDG.h"
  51 #include "TParticle.h"
  52 #include "TString.h"
  53
  54 ClassImp(AliMUONTrackLight)
  55
  56 //===================================================================
  57
  58 AliMUONTrackLight::AliMUONTrackLight()
  59   : TObject(),
  60     fPrec(),
  61     fIsTriggered(kFALSE),
  62     fCharge(-999),
  63     fChi2(-1),
  64     fCentr(-1),
  65     fPgen(),
  66     fTrackPythiaLine(-999),
  67     fTrackPDGCode(-999),
  68     fOscillation(kFALSE),
  69     fNParents(0),
  70     fWeight(1)
  71 {
  72   /// default constructor
  73   fPgen.SetPxPyPzE(0.,0.,0.,0.);
  74   fPrec.SetPxPyPzE(0.,0.,0.,0.);
  75   for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
  76   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
  77     fParentPDGCode[npar] = -1;
  78     fParentPythiaLine[npar] = -1;
  79   }
  80   for (Int_t i = 0; i < 4; i++){
  81     fQuarkPDGCode[i] = -1;
  82     fQuarkPythiaLine[i] = -1;
  83   }
  84 }
  85
  86 //============================================
  87 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
  88   : TObject(muonCopy),
  89     fPrec(muonCopy.fPrec),
  90     fIsTriggered(muonCopy.fIsTriggered),
  91     fCharge(muonCopy.fCharge),
  92     fChi2(muonCopy.fChi2),
  93     fCentr(muonCopy.fCentr),
  94     fPgen(muonCopy.fPgen),
  95     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
  96     fTrackPDGCode(muonCopy.fTrackPDGCode),
  97     fOscillation(muonCopy.fOscillation),
  98     fNParents(muonCopy.fNParents),
  99     fWeight(muonCopy.fWeight)
 100 {
 101   /// copy constructor
 102   for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
 103   for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
 104     fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
 105     fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
 106   }
 107   for (Int_t i = 0; i < 4; i++){
 108     fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
 109     fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
 110   }
 111 }
 112
 113 //============================================
 114 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
 115   : TObject(),
 116     fPrec(),
 117     fIsTriggered(kFALSE),
 118     fCharge(-999),
 119     fChi2(-1),
 120     fCentr(-1),
 121     fPgen(),
 122     fTrackPythiaLine(-999),
 123     fTrackPDGCode(-999),
 124     fOscillation(kFALSE),
 125     fNParents(0),
 126     fWeight(1)
 127 {
 128   /// constructor
 129   //AliMUONTrackLight();
 130   FillFromESD(muonTrack);
 131 }
 132
 133 //============================================
 134 AliMUONTrackLight::~AliMUONTrackLight()
 135 {
 136 /// Destructor
 137 }
 138
 139 //============================================
 140
 141 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
 142   /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
 143   AliMUONTrackParam trPar(*((AliMUONTrackParam*) (trackReco->GetTrackParamAtCluster()->First())));
 144   //  AliMUONTrackParam *trPar  = trackReco->GetTrackParamAtVertex();
 145   AliMUONTrackExtrap::ExtrapToVertex(&trPar,0.,0.,0.,0.,0.);
 146   this->SetCharge(Int_t(TMath::Sign(1.,trPar.GetInverseBendingMomentum())));
 147   this->SetPxPyPz(trPar.Px(),trPar.Py(), trPar.Pz());
 148   this->SetTriggered(trackReco->GetMatchTrigger());
 149
 150   Double_t xyz[3] = { trPar.GetNonBendingCoor(),
 151                       trPar.GetBendingCoor(),
 152                       zvert};
 153   this->SetVertex(xyz);
 154 }
 155
 156 //============================================
 157 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
 158   /// computes prec and charge from ESD track
 159   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 160   Double_t thetaX = muonTrack->GetThetaX();
 161   Double_t thetaY = muonTrack->GetThetaY();
 162   Double_t tanthx = TMath::Tan(thetaX);
 163   Double_t tanthy = TMath::Tan(thetaY);
 164   Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
 165   Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
 166   Double_t px     = pz * tanthx;
 167   Double_t py     = pz * tanthy;
 168   fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
 169   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 170   fPrec.SetPxPyPzE(px,py,pz,energy);
 171   // get the position
 172   fXYZ[0] = muonTrack->GetNonBendingCoor();
 173   fXYZ[1] = muonTrack->GetBendingCoor();
 174   if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
 175   else fXYZ[2] = zvert;
 176   // get the chi2 per d.o.f.
 177   fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
 178   fIsTriggered = muonTrack->GetMatchTrigger();
 179 }
 180
 181 //============================================
 182 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
 183   /// set the reconstructed 4-momentum, assuming the particle is a muon
 184   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 185   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
 186   fPrec.SetPxPyPzE(px,py,pz,energy);
 187 }
 188
 189 //============================================
 190 TParticle* AliMUONTrackLight::FindRefTrack(
 191                 AliMUONTrack* trackReco, AliMUONVTrackStore* trackRefArray,
 192                 AliStack* stack, Bool_t tempfix
 193         )
 194 {
 195   /// Find the Monte Carlo (MC) particle that corresponds to a given reconstructed track.
 196   /// @param trackReco  This is the reconstructed track for which we want to find a
 197   ///           corresponding MC particle.
 198   /// @param trackRefArray  The list of MC reference tracks as generated by
 199   ///           AliMUONRecoCheck::ReconstructedTracks for example.
 200   /// @param stack  The MC stack of simulated particles.
 201   /// @param tempfix  This specifies if a temporary fix should be applied, where we only
 202   ///           check if the first hit on chamber 1, the vertex momentum, the X-Y vertex
 203   ///           coordinate and X-Y slope match between 'trackReco' and a MC reference
 204   ///           track in 'trackRefArray'.
 205   ///           This temporary fix is necessary while AliMUONTrack objects are no longer
 206   ///           stored after reconstruction and the hit information is not completely
 207   ///           stored in the ESD data.
 208
 209   TParticle *part = 0;
 210   const Double_t kSigma2Cut = 16;  // 4 sigmas cut, kSigma2Cut = 4*4
 211   Int_t compPart = 0;
 212   TIter next(trackRefArray->CreateIterator());
 213   AliMUONTrack* trackRef;
 214   while ( (trackRef = static_cast<AliMUONTrack*>(next())) ) {
 215     // check if trackRef is compatible with trackReco:
 216     //routine returns for each chamber a yes/no information if the
 217     //hit of rec. track and hit of referenced track are compatible
 218     Bool_t *compTrack = trackRef->CompatibleTrack(trackReco,kSigma2Cut);
 219     Int_t iTrack = 0;
 220     if (tempfix)
 221     {
 222       // All we can check with this temporary fix is that the first hit on chamber 1
 223       // is the same and that the particles momentum, slope and vertex X and Y are
 224       // compatible within resolution.
 225       if (compTrack[0]) iTrack++;
 226       AliMUONTrackParam* paramA = trackRef->GetTrackParamAtVertex();
 227       AliMUONTrackParam* paramB = trackReco->GetTrackParamAtVertex();
 228       if (paramA != NULL and paramB != NULL)
 229       {
 230         // Fetch the variances.
 231         Double_t varX = paramA->GetCovariances()[0][0] + paramB->GetCovariances()[0][0];
 232         Double_t varSlopeX = paramA->GetCovariances()[1][1] + paramB->GetCovariances()[1][1];
 233         Double_t varY = paramA->GetCovariances()[2][2] + paramB->GetCovariances()[2][2];
 234         Double_t varSlopeY = paramA->GetCovariances()[3][3] + paramB->GetCovariances()[3][3];
 235         Double_t varInvP = paramA->GetCovariances()[4][4] + paramB->GetCovariances()[4][4];
 236         // We might have to fill these variances with defaults.
 237         if (varX == 0) varX = AliMUONConstants::DefaultNonBendingReso2();
 238         if (varSlopeX == 0) varSlopeX = AliMUONConstants::DefaultNonBendingReso2();
 239         if (varY == 0) varY = AliMUONConstants::DefaultBendingReso2();
 240         if (varSlopeY == 0) varSlopeY = AliMUONConstants::DefaultBendingReso2();
 241         if (varInvP == 0)
 242         {
 243                 // Use ~10% resolution for momentum.
 244                 Double_t sigmaA = 0.1 * paramA->GetInverseBendingMomentum();
 245                 Double_t sigmaB = 0.1 * paramB->GetInverseBendingMomentum();
 246                 varInvP = sigmaA*sigmaA + sigmaB*sigmaB;
 247         }
 248         // Calculate the differences between parameters.
 249         Double_t diffX = paramA->GetNonBendingCoor() - paramB->GetNonBendingCoor();
 250         Double_t diffSlopeX = paramA->GetNonBendingSlope() - paramB->GetNonBendingSlope();
 251         Double_t diffY = paramA->GetBendingCoor() - paramB->GetBendingCoor();
 252         Double_t diffSlopeY = paramA->GetBendingSlope() - paramB->GetBendingSlope();
 253         Double_t diffInvP = paramA->GetInverseBendingMomentum() - paramB->GetInverseBendingMomentum();
 254         // Check that all the parameters are within kSigma2Cut limit.
 255         if (diffX*diffX / varX < kSigma2Cut and diffY*diffY / varY < kSigma2Cut) iTrack++;
 256         if (diffSlopeX*diffSlopeX / varSlopeX < kSigma2Cut and diffSlopeY*diffSlopeY / varSlopeY < kSigma2Cut) iTrack++;
 257         if (diffInvP*diffInvP / varInvP < kSigma2Cut) iTrack++;
 258       }
 259     }
 260     else
 261       iTrack = this->TrackCheck(compTrack); //returns number of validated conditions
 262
 263     if (iTrack==4) {
 264       compPart++;
 265       Int_t trackID = trackRef->GetTrackID();
 266       this->SetTrackPythiaLine(trackID);
 267       part = stack->Particle(trackID);
 268       fTrackPDGCode = part->GetPdgCode();
 269     }
 270   }
 271   if (compPart>1) {
 272     AliFatal("More than one particle compatible with the reconstructed track.");
 273   }
 274   return part;
 275 }
 276
 277 //============================================
 278 Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
 279   /// Apply reconstruction requirements
 280   /// Return number of validated conditions
 281   /// If all the tests are verified then TrackCheck = 4 (good track)
 282
 283   Int_t iTrack = 0;
 284   Int_t hitsInLastStations = 0;
 285
 286   // apply reconstruction requirements
 287   if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
 288   if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
 289   if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
 290   for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
 291     if (compTrack[ch]) hitsInLastStations++;
 292   }
 293   if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
 294   return iTrack;
 295 }
 296
 297 //============================================
 298 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
 299   /// scans the muon history to determine parents pdg code and pythia line
 300   Int_t countP = -1;
 301   Int_t parents[10], parLine[10];
 302   Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
 303
 304   TParticle *mother;
 305   Int_t status=-1, pdg=-1;
 306   while(lineM >= 0){
 307
 308     mother = stack->Particle(lineM); //direct mother of rec. track
 309     pdg = mother->GetPdgCode();//store PDG code of first mother
 310     // break if a string, gluon, quark or diquark is found
 311     if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
 312     parents[++countP] = pdg;
 313     parLine[countP] = lineM;
 314     status = mother->GetStatusCode();//get its status code to check if oscillation occured
 315     if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
 316     lineM = mother->GetFirstMother();
 317   }
 318   //store all the fragmented parents in an array:
 319   for(int i = 0; i <= countP; i++){
 320     this->SetParentPDGCode(i,parents[countP-i]);
 321     this->SetParentPythiaLine(i,parLine[countP-i]);
 322   }
 323   fNParents = countP+1;
 324   countP = -1;
 325
 326   //and store the lines of the string and further quarks in another array:
 327   while(lineM >= 0){
 328     mother = stack->Particle(lineM);
 329     pdg = mother->GetPdgCode();
 330     //now, get information before the fragmentation
 331     this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
 332     this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
 333     lineM = mother->GetFirstMother();
 334   }
 335
 336   //check if in case of HF production, the string points to the correct end
 337   //and correct it in case of need:
 338   countP = 1;
 339   for(int par = 0; par < 4; par++){
 340     if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
 341       countP = par; //get the quark just before hadronisation
 342       break;
 343     }
 344   }
 345   if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
 346     if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
 347
 348       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 349           this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
 350         );
 351
 352       Int_t pdg = this->GetQuarkPDGCode(countP), line = this->GetQuarkPythiaLine(countP);
 353       this->ResetQuarkInfo();
 354       while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
 355                                                         //must coincide with the flavour of the last fragmented mother
 356
 357         pdg = stack->Particle(++line)->GetPdgCode();
 358       }
 359       //now, we have the correct string end and the correct line
 360       //continue to fill again all parents and corresponding lines
 361       while(line >= 0){
 362         mother = stack->Particle(line);//get again the mother
 363         pdg = mother->GetPdgCode();
 364         this->SetQuarkPythiaLine(countP, line);
 365         this->SetQuarkPDGCode(countP++, pdg);
 366         line = mother->GetFirstMother();
 367       }
 368       this->PrintInfo("h");
 369     }//mismatch
 370   }
 371 }
 372
 373 //====================================
 374 void AliMUONTrackLight::ResetQuarkInfo(){
 375   /// resets parton information
 376   for(int pos = 1; pos < 4; pos++){//[0] is the string
 377     this->SetQuarkPDGCode(pos,-1);
 378     this->SetQuarkPythiaLine(pos,-1);
 379   }
 380 }
 381
 382 //====================================
 383 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
 384   /// checks if the particle is a B0
 385   Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
 386   Bool_t answer = kFALSE;
 387   for(int i = 0; i < 2; i++){
 388     if(TMath::Abs(intTest) == bMes0[i]){
 389       answer = kTRUE;
 390       break;
 391     }
 392   }
 393   return answer;
 394 }
 395 //====================================
 396 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
 397   /// checks if mother is a resonance
 398   Int_t intTest = GetParentPDGCode(index);
 399   // the resonance pdg code is built this way
 400   // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
 401   Int_t id=intTest%100000;
 402   return (!((id-id%10)%110));
 403 }
 404 //====================================
 405 Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
 406   /// returns the flavour of parent idParent (idParent=0 is the oldest
 407   /// hadronized parent)
 408   Int_t pdg = GetParentPDGCode(idParent);
 409   Int_t quark = TMath::Abs(pdg/100);
 410   if(quark > 9) quark = quark/10;
 411   return quark;
 412 }
 413
 414 //====================================
 415 void AliMUONTrackLight::PrintInfo(Option_t* opt){
 416   /// prints information about the track:
 417   /// - "H" muon's decay history
 418   /// - "K" muon kinematics
 419   /// - "A" all variables
 420   TString options(opt);
 421   options.ToUpper();
 422
 423   if(options.Contains("H") || options.Contains("A")){ //muon decay history
 424     char *name= new char[100];
 425     TString pdg = "", line = "";
 426     for(int i = 3; i >= 0; i--){
 427       if(this->GetQuarkPythiaLine(i)>= 0){
 428         sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
 429         line += name;
 430         sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
 431         pdg += name;
 432       }
 433     }
 434     for(int i = 0; i < fNParents; i++){
 435       if(this->GetParentPythiaLine(i)>= 0){
 436         sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
 437         line += name;
 438         sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
 439         pdg += name;
 440       }
 441     }
 442     sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
 443     sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
 444
 445     printf("\nmuon's decay history:\n");
 446     printf(" PDG: %s\n", pdg.Data());
 447     printf("line: %s\n", line.Data());
 448   }
 449   if(options.Contains("K") || options.Contains("A")){ //muon kinematic
 450
 451     Int_t charge = this->GetCharge();
 452     Double_t *vtx = this->GetVertex();
 453     TLorentzVector momRec = this->GetPRec();
 454     TLorentzVector momGen = this->GetPGen();
 455     printf("the track's charge is %d\n", charge);
 456     printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 457     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 458     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 459     printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
 460            momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 461            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 462   }
 463 }
 464 //====================================
 465 Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
 466   /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
 467   Int_t pdg = this->GetParentPDGCode(idparent);
 468   if (TMath::Abs(pdg)==211 || //pi+
 469       TMath::Abs(pdg)==321 || //K+
 470       TMath::Abs(pdg)==213 || //rho+
 471       TMath::Abs(pdg)==311 || //K0
 472       TMath::Abs(pdg)==313 || //K*0
 473       TMath::Abs(pdg)==323    //K*+
 474       ) {
 475     return kTRUE;
 476   }
 477   else return kFALSE;
 478 }
 479 //====================================
 480 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg){
 481   /// check if the provided pdg code corresponds to a diquark
 482   pdg = TMath::Abs(pdg);
 483   if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
 484   else return kFALSE;
 485 }