PWG3/muon/AliMuonInfoStoreMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   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 purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // class used to extract and store info of MC particle
  19 //
  20 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
  21 //                     zhangxm@iopp.ccnu.edu.cn
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25
  26 #include "AliMCEvent.h"
  27 #include "AliAODMCParticle.h"
  28 #include "AliESDMuonTrack.h"
  29 #include "AliAODTrack.h"
  30 #include "AliMuonInfoStoreRD.h"
  31 #include "AliMuonInfoStoreMC.h"
  32 #include "AliMCEvent.h"
  33
  34 ClassImp(AliMuonInfoStoreMC)
  35
  36 const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
  37 const Int_t   AliMuonInfoStoreMC::fgkSourcesN = 6;
  38
  39 //-----------------------------------------------------------------------------
  40 AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
  41 AliMuonInfoStoreRD(),
  42 fIsFull(kFALSE),
  43 fLorentzP(),
  44 fTrackIndex(-1),
  45 fTrackPDGCode(0),
  46 fSource(-1),
  47 fNParents(0),
  48 fOscillation(kFALSE),
  49 fWeight(0.)
  50 {
  51   //
  52   // default constructor
  53   //
  54   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  55   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  56 }
  57
  58 //-----------------------------------------------------------------------------
  59 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, AliMCEvent *mcEvent, Bool_t full) :
  60 AliMuonInfoStoreRD(trkAOD),
  61 fIsFull(full),
  62 fLorentzP(),
  63 fTrackIndex(-1),
  64 fTrackPDGCode(0),
  65 fSource(-1),
  66 fNParents(0),
  67 fOscillation(kFALSE),
  68 fWeight(0.)
  69 {
  70   //
  71   // default constructor
  72   //
  73   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  74   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  75
  76   this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
  77 }
  78
  79 //-----------------------------------------------------------------------------
  80 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEvent *mcEvent, Bool_t full) :
  81 AliMuonInfoStoreRD(trkESD),
  82 fIsFull(full),
  83 fLorentzP(),
  84 fTrackIndex(-1),
  85 fTrackPDGCode(0),
  86 fSource(-1),
  87 fNParents(0),
  88 fOscillation(kFALSE),
  89 fWeight(0.)
  90 {
  91   //
  92   // default constructor
  93   //
  94   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  95   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  96
  97   this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
  98 }
  99
 100 //-----------------------------------------------------------------------------
 101 AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
 102 AliMuonInfoStoreRD(src),
 103 fIsFull(src.fIsFull),
 104 fLorentzP(src.fLorentzP),
 105 fTrackIndex(src.fTrackIndex),
 106 fTrackPDGCode(src.fTrackPDGCode),
 107 fSource(src.fSource),
 108 fNParents(src.fNParents),
 109 fOscillation(src.fOscillation),
 110 fWeight(src.fWeight)
 111 {
 112   //
 113   // copy constructor
 114   //
 115   for (Int_t i=5; i--;) {
 116     fParentIndex[i]   = src.fParentIndex[i];
 117     fParentPDGCode[i] = src.fParentPDGCode[i];
 118   }
 119   for (Int_t i=4; i--;) {
 120     fQuarkIndex[i]    = src.fQuarkIndex[i];
 121     fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
 122   }
 123 }
 124
 125 //-----------------------------------------------------------------------------
 126 AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
 127 {
 128   //
 129   // assignment constructor
 130   //
 131   if(&src==this) return *this;
 132   AliMuonInfoStoreRD::operator=(src);
 133
 134   fIsFull       = src.fIsFull;
 135   fLorentzP     = src.fLorentzP;
 136   fTrackIndex   = src.fTrackIndex;
 137   fTrackPDGCode = src.fTrackPDGCode;
 138   fSource       = src.fSource;
 139   fNParents     = src.fNParents;
 140   fOscillation  = src.fOscillation;
 141   fWeight       = src.fWeight;
 142   for (Int_t i=5; i--;) {
 143     fParentIndex[i]   = src.fParentIndex[i];
 144     fParentPDGCode[i] = src.fParentPDGCode[i];
 145   }
 146   for (Int_t i=4; i--;) {
 147     fQuarkIndex[i]    = src.fQuarkIndex[i];
 148     fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
 149   }
 150
 151   return *this;
 152 }
 153
 154 //-----------------------------------------------------------------------------
 155 AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
 156 {
 157   //
 158   // destructor
 159   //
 160 }
 161
 162 //-----------------------------------------------------------------------------
 163 void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
 164 {
 165   // fill track MC info with AOD base
 166   fTrackIndex = label;
 167   if (fTrackIndex<0) { fSource=5; return; }
 168
 169   AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
 170   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
 171
 172   fTrackPDGCode = pMC->PdgCode();
 173   if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
 174
 175   Int_t lineM = pMC->GetMother();
 176   if (lineM<0) { fSource=2; return; }
 177
 178   Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
 179   if (!isPrimary) { fSource=3; return; }
 180
 181   Int_t countP=-1, pdg=0;
 182   Int_t parents[10], parLine[10];
 183   AliAODMCParticle *mother = 0;
 184   while(lineM>=0){
 185     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
 186     pdg = mother->GetPdgCode();
 187     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 188     parents[++countP] = pdg;
 189     parLine[countP] = lineM;
 190     lineM = mother->GetMother();
 191   }
 192   for(Int_t i=0; i<=countP; i++){
 193     fParentIndex[i] = parLine[countP-i];
 194     fParentPDGCode[i] = parents[countP-i];
 195   }
 196   fNParents = countP + 1;
 197
 198   if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
 199
 200   fSource = this->SelectHFMuon();
 201   return;
 202 }
 203
 204 //-----------------------------------------------------------------------------
 205 void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
 206 {
 207   // fill track MC info with ESD base
 208   fTrackIndex = label;
 209   if (fTrackIndex<0) { fSource=5; return; }
 210
 211   TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
 212   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
 213
 214   fTrackPDGCode = pMC->GetPdgCode();
 215   if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
 216
 217   Int_t lineM = pMC->GetFirstMother();
 218   if (lineM<0) { fSource=2; return; }
 219
 220   if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
 221
 222   Int_t countP=-1, pdg=0;
 223   Int_t parents[10], parLine[10];
 224   TParticle *mother = 0;
 225   while(lineM>=0){
 226     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
 227     pdg = mother->GetPdgCode();
 228     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 229     parents[++countP] = pdg;
 230     parLine[countP] = lineM;
 231     lineM = mother->GetFirstMother();
 232   }
 233   for(Int_t i=0; i<=countP; i++){
 234     fParentIndex[i] = parLine[countP-i];
 235     fParentPDGCode[i] = parents[countP-i];
 236   }
 237   fNParents = countP + 1;
 238
 239   if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
 240
 241   fSource = this->SelectHFMuon();
 242   return;
 243 }
 244
 245 //-----------------------------------------------------------------------------
 246 void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
 247 {
 248   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 249
 250   if (lineM<0) return;
 251   Int_t countP=-1, pdg=0;
 252   AliAODMCParticle *mother = 0;
 253   while(lineM>=0){
 254     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
 255     pdg = mother->GetPdgCode();
 256     fQuarkIndex[++countP] = lineM;
 257     fQuarkPDGCode[countP] = pdg;
 258     lineM = mother->GetMother();
 259   }
 260
 261   // for PYTHIA checking
 262   countP = 1;
 263   for(Int_t par=0; par<4; par++) {
 264     if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
 265   }
 266   if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
 267     if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
 268       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 269                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 270
 271       pdg = this->QuarkPDGCode(countP);
 272       Int_t line = this->QuarkIndex(countP);
 273       this->ResetQuarkInfo();
 274       while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
 275         pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
 276       }
 277       while(line>=0){
 278         mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
 279         pdg = mother->GetPdgCode();
 280         fQuarkIndex[countP] = line;
 281         fQuarkPDGCode[countP++] = pdg;
 282         line = mother->GetMother();
 283       }
 284     }
 285   }
 286   return;
 287 }
 288
 289 //-----------------------------------------------------------------------------
 290 void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
 291 {
 292   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 293
 294   if (lineM<0) return;
 295   Int_t countP=-1, pdg=0;
 296   TParticle *mother = 0;
 297   while(lineM>=0){
 298     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
 299     pdg = mother->GetPdgCode();
 300     fQuarkIndex[++countP] = lineM;
 301     fQuarkPDGCode[countP] = pdg;
 302     lineM = mother->GetFirstMother();
 303   }
 304
 305   // for PYTHIA checking
 306   countP = 1;
 307   for(Int_t par=0; par<4; par++) {
 308     if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
 309   }
 310   if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
 311     if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
 312       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 313                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 314
 315       pdg = this->QuarkPDGCode(countP);
 316       Int_t line = this->QuarkIndex(countP);
 317       this->ResetQuarkInfo();
 318       while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
 319         pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
 320       }
 321       while(line>=0){
 322         mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
 323         pdg = mother->GetPdgCode();
 324         fQuarkIndex[countP] = line;
 325         fQuarkPDGCode[countP++] = pdg;
 326         line = mother->GetFirstMother();
 327       }
 328     }
 329   }
 330   return;
 331 }
 332
 333 //-----------------------------------------------------------------------------
 334 Int_t AliMuonInfoStoreMC::SelectHFMuon()
 335 {
 336   // set info of muon from HF
 337
 338   Int_t flv = ParentFlavour(0);
 339   if (flv!=4 && flv!=5) return 2;
 340
 341   Bool_t isRes = kFALSE;
 342   Int_t i=0, nparents=this->ParentsN();
 343   while (i<nparents && !isRes) {
 344     isRes = IsMotherAResonance(i++);
 345   }
 346
 347   if (isRes) return 2;
 348   if (flv==5) return 0;
 349   else return 1;
 350 }
 351
 352 //-----------------------------------------------------------------------------
 353 Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
 354 {
 355   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 356   pdg = TMath::Abs(pdg);
 357   if(pdg>1000 && (pdg%100)<10) return kTRUE;
 358   else return kFALSE;
 359 }
 360
 361 //-----------------------------------------------------------------------------
 362 void AliMuonInfoStoreMC::ResetQuarkInfo()
 363 {
 364   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 365   for(Int_t pos=1; pos<4; pos++) {
 366     fQuarkIndex[pos] = -1;
 367     fQuarkPDGCode[pos] = 0;
 368   }
 369   return;
 370 }
 371
 372 //-----------------------------------------------------------------------------
 373 Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
 374 {
 375   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 376   Int_t pdg = ParentPDGCode(i);
 377   pdg = TMath::Abs(pdg/100);
 378   if(pdg>9) pdg /= 10;
 379   return pdg;
 380 }
 381
 382 //-----------------------------------------------------------------------------
 383 Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
 384 {
 385   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 386   Int_t pdg = ParentPDGCode(i);
 387   Int_t id=pdg%100000;
 388   return (!((id-id%10)%110));
 389 }