PWG3/dielectron/AliDielectronMC.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 AliDielectronMC                #
  19 //#       Cut Class for Jpsi->e+e- analysis           #
  20 //#                                                   #
  21 //#   by WooJin J. Park, GSI / W.J.Park@gsi.de        #
  22 //#                                                   #
  23 //#####################################################
  24
  25 #include <AliAnalysisManager.h>
  26 #include <AliAODHandler.h>
  27 #include <AliESDInputHandler.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliMCParticle.h>
  31 #include <AliAODMCParticle.h>
  32 #include <AliStack.h>
  33 #include <AliESDEvent.h>
  34 #include <AliESDtrack.h>
  35 #include <AliLog.h>
  36
  37 #include "AliDielectronMC.h"
  38
  39 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
  40
  41 //____________________________________________________________
  42 AliDielectronMC* AliDielectronMC::Instance()
  43 {
  44   //
  45   // return pointer to singleton implementation
  46   //
  47   if (fgInstance) return fgInstance;
  48
  49   AnalysisType type=kUNSET;
  50   Bool_t hasMC=kFALSE;
  51   if (AliAnalysisManager::GetAnalysisManager()){
  52     if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
  53     else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD;
  54
  55     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  56     hasMC=mcHandler!=0x0;
  57   }
  58
  59   fgInstance=new AliDielectronMC(type);
  60   fgInstance->SetHasMC(hasMC);
  61
  62   return fgInstance;
  63 }
  64
  65 //____________________________________________________________
  66 AliDielectronMC::AliDielectronMC(AnalysisType type):
  67   fMCEvent(0x0),
  68   fStack(0x0),
  69   fAnaType(type),
  70   fHasMC(kTRUE)
  71 {
  72   //
  73   // default constructor
  74   //
  75 }
  76
  77
  78 //____________________________________________________________
  79 AliDielectronMC::~AliDielectronMC()
  80 {
  81   //
  82   // default destructor
  83   //
  84
  85 }
  86
  87 //____________________________________________________________
  88 void AliDielectronMC::Initialize()
  89 {
  90   //
  91   // initialize MC class
  92   //
  93   if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
  94 }
  95
  96 //____________________________________________________________
  97 Int_t AliDielectronMC::GetNMCTracks()
  98 {
  99   //
 100   //  return the number of generated tracks from MC event
 101   //
 102   if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
 103   return fMCEvent->GetNumberOfTracks();
 104 }
 105
 106 //____________________________________________________________
 107 Int_t AliDielectronMC::GetNMCTracksFromStack()
 108 {
 109   //
 110   //  return the number of generated tracks from stack
 111   //
 112   if (!fStack){ AliError("No fStack"); return -999; }
 113   return fStack->GetNtrack();
 114 }
 115
 116 //____________________________________________________________
 117 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk)
 118 {
 119   //
 120   // return MC track directly from MC event
 121   //
 122   if (itrk<0) return NULL;
 123   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 124   AliVParticle * track = fMCEvent->GetTrack(itrk); //  tracks from MC event
 125   return track;
 126 }
 127
 128 //____________________________________________________________
 129 Bool_t AliDielectronMC::ConnectMCEvent()
 130 {
 131   //
 132   // connect stack object from the mc handler
 133   //
 134   fMCEvent=0x0;
 135   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 136   if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
 137   if (!mcHandler->InitOk() ) return kFALSE;
 138   if (!mcHandler->TreeK() )  return kFALSE;
 139   if (!mcHandler->TreeTR() ) return kFALSE;
 140
 141   AliMCEvent* mcEvent = mcHandler->MCEvent();
 142   if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
 143   fMCEvent = mcEvent;
 144
 145   if (!UpdateStack()) return kFALSE;
 146   return kTRUE;
 147 }
 148
 149 //____________________________________________________________
 150 Bool_t AliDielectronMC::UpdateStack()
 151 {
 152   //
 153   // update stack with new event
 154   //
 155   if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
 156   AliStack* stack = fMCEvent->Stack();
 157   if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
 158   fStack = stack;
 159   return kTRUE;
 160 }
 161
 162 //____________________________________________________________
 163 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
 164 {
 165   //
 166   // return MC track
 167   //
 168   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 169
 170   Int_t nStack = fMCEvent->GetNumberOfTracks();
 171   Int_t label = TMath::Abs(_track->GetLabel());
 172   if(label>nStack)return NULL;
 173
 174   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
 175   return mctrack;
 176 }
 177
 178 //____________________________________________________________
 179 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
 180 {
 181   //
 182   // return MC track from stack
 183   //
 184   Int_t label = TMath::Abs(_track->GetLabel());
 185   if (!fStack) AliWarning("fStack is not available. Update stack first.");
 186   TParticle* mcpart = fStack->Particle(label);
 187   if (!mcpart) return NULL;
 188   return mcpart;
 189 }
 190
 191 //____________________________________________________________
 192 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
 193 {
 194   //
 195   // return MC track mother
 196   //
 197   AliMCParticle* mcpart = GetMCTrack(_track);
 198   if (!mcpart) return NULL;
 199   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
 200   if (!mcmother) return NULL;
 201   return mcmother;
 202 }
 203
 204 //____________________________________________________________
 205 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
 206   //
 207   // return MC track mother
 208   //
 209   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
 210   return mcmother;
 211 }
 212
 213 //____________________________________________________________
 214 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
 215   //
 216   // return MC track mother
 217   //
 218   AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
 219   return mcmother;
 220 }
 221
 222 //____________________________________________________________
 223 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
 224 {
 225   //
 226   // return MC track mother from stack
 227   //
 228   TParticle* mcpart = GetMCTrackFromStack(_track);
 229   if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL;
 230   TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
 231   if (!mcmother) return NULL;
 232   return mcmother;
 233 }
 234
 235 //____________________________________________________________
 236 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
 237 {
 238   //
 239   // return PDG code of the track from the MC truth info
 240   //
 241   AliMCParticle* mcpart = GetMCTrack(_track);
 242   if (!mcpart) return -999;
 243   return mcpart->PdgCode();
 244 }
 245
 246 //____________________________________________________________
 247 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
 248 {
 249   //
 250   // return MC PDG code from stack
 251   //
 252   TParticle* mcpart = GetMCTrackFromStack(_track);
 253   if (!mcpart) return -999;
 254   return mcpart->GetPdgCode();
 255 }
 256
 257 //____________________________________________________________
 258 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
 259 {
 260   //
 261   // return PDG code of the mother track from the MC truth info
 262   //
 263   AliMCParticle* mcmother = GetMCTrackMother(_track);
 264   if (!mcmother) return -999;
 265   return mcmother->PdgCode();
 266 }
 267
 268 //____________________________________________________________
 269 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
 270 {
 271   //
 272   // return PDG code of the mother track from stack
 273   //
 274   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 275   if (!mcmother) return -999;
 276   return mcmother->GetPdgCode();
 277 }
 278
 279 //____________________________________________________________
 280 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
 281 {
 282   //
 283   // return process number of the track
 284   //
 285   AliMCParticle* mcpart = GetMCTrack(_track);
 286   if (!mcpart) return -999;
 287   return 0;
 288 }
 289
 290 //____________________________________________________________
 291 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
 292 {
 293   //
 294   // return process number of the track
 295   //
 296   TParticle* mcpart = GetMCTrackFromStack(_track);
 297   if (!mcpart) return -999;
 298   return mcpart->GetUniqueID();
 299 }
 300
 301 //____________________________________________________________
 302 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
 303 {
 304   //
 305   // returns the number of daughters
 306   //
 307   AliMCParticle *mcmother=GetMCTrackMother(track);
 308   if(!mcmother||!mcmother->Particle()) return -999;
 309 //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 310   return mcmother->Particle()->GetNDaughters();
 311 }
 312
 313 //____________________________________________________________
 314 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
 315 {
 316   //
 317   // returns the number of daughters
 318   //
 319   AliMCParticle *mcmother=GetMCTrackMother(particle);
 320   if(!mcmother||!mcmother->Particle()) return -999;
 321 //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 322   return mcmother->Particle()->GetNDaughters();
 323 }
 324
 325 //____________________________________________________________
 326 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
 327 {
 328   //
 329   // returns the number of daughters
 330   //
 331   AliAODMCParticle *mcmother=GetMCTrackMother(particle);
 332   if(!mcmother) return -999;
 333   return mcmother->GetNDaughters();
 334 }
 335
 336 //____________________________________________________________
 337 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
 338 {
 339   //
 340   // return process number of the mother of the track
 341   //
 342   AliMCParticle* mcmother = GetMCTrackMother(_track);
 343   if (!mcmother) return -999;
 344   return 0;
 345 }
 346
 347 //____________________________________________________________
 348 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
 349 {
 350   //
 351   // return process number of the mother of the track
 352   //
 353   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 354   if (!mcmother) return -999;
 355   return mcmother->GetUniqueID();
 356 }
 357
 358 //____________________________________________________________
 359 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
 360 {
 361   //
 362   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 363   //
 364
 365   if (!fMCEvent) return kFALSE;
 366   if (!particle) return kFALSE;
 367
 368   if (particle->IsA()==AliMCParticle::Class()){
 369     return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
 370   } else if (particle->IsA()==AliAODMCParticle::Class()){
 371     return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
 372   } else {
 373     AliError("Unknown particle type");
 374   }
 375   return kFALSE;
 376
 377 }
 378
 379 //____________________________________________________________
 380 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
 381 {
 382   //
 383   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 384   // ESD case
 385   //
 386
 387   //check pdg code
 388   if (particle->PdgCode()!=pdgMother) return kFALSE;
 389   Int_t ifirst = particle->GetFirstDaughter();
 390   Int_t ilast  = particle->GetLastDaughter();
 391
 392   //check number of daughters
 393   if ((ilast-ifirst)!=1) return kFALSE;
 394   AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 395   AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
 396
 397   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 398   if (firstD->Charge()>0){
 399     if (firstD->PdgCode()!=-11) return kFALSE;
 400     if (secondD->PdgCode()!=11) return kFALSE;
 401   }else{
 402     if (firstD->PdgCode()!=11) return kFALSE;
 403     if (secondD->PdgCode()!=-11) return kFALSE;
 404   }
 405
 406   return kTRUE;
 407 }
 408
 409 //____________________________________________________________
 410 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
 411 {
 412   //
 413   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 414   // AOD case
 415   //
 416   if (particle->GetPdgCode()!=pdgMother) return kFALSE;
 417   if (particle->GetNDaughters()!=2) return kFALSE;
 418
 419   Int_t ifirst = particle->GetDaughter(0);
 420   Int_t ilast  = particle->GetDaughter(1);
 421
 422   //check number of daughters
 423   if ((ilast-ifirst)!=1) return kFALSE;
 424   AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 425   AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
 426
 427   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 428   if (firstD->Charge()>0){
 429     if (firstD->GetPdgCode()!=11) return kFALSE;
 430     if (secondD->GetPdgCode()!=-11) return kFALSE;
 431   }else{
 432     if (firstD->GetPdgCode()!=-11) return kFALSE;
 433     if (secondD->GetPdgCode()!=11) return kFALSE;
 434   }
 435   return kTRUE;
 436 }
 437
 438 //____________________________________________________________
 439 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 440 {
 441   //
 442   // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
 443   //
 444   if (!fMCEvent) return -1;
 445
 446   if (fAnaType==kESD) return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
 447   else if (fAnaType==kAOD) return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
 448
 449   return -1;
 450 }
 451
 452 //____________________________________________________________
 453 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 454 {
 455   //
 456   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 457   //    have the same mother and the mother had pdg code pdgMother
 458   // ESD case
 459   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 460   //
 461
 462   AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
 463   AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
 464
 465   if (!mcPart1||!mcPart2) return -1;
 466
 467   Int_t lblMother1=mcPart1->GetMother();
 468   Int_t lblMother2=mcPart2->GetMother();
 469
 470   AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 471   if (!mcMother1) return -1;
 472   if (lblMother1!=lblMother2) return -1;
 473   if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
 474   if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
 475   if (mcMother1->PdgCode()!=pdgMother) return -1;
 476
 477   return lblMother1;
 478 }
 479
 480 //____________________________________________________________
 481 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 482 {
 483   //
 484   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 485   //    have the same mother and the mother had pdg code pdgMother
 486   // AOD case
 487   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 488   //
 489   AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
 490   AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
 491
 492   if (!mcPart1||!mcPart2) return -1;
 493
 494   Int_t lblMother1=mcPart1->GetMother();
 495   Int_t lblMother2=mcPart2->GetMother();
 496
 497   AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 498
 499   if (!mcMother1) return -1;
 500   if (lblMother1!=lblMother2) return -1;
 501   if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
 502   if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
 503   if (mcMother1->GetPdgCode()!=pdgMother) return -1;
 504
 505   return lblMother1;
 506 }
 507
 508 //____________________________________________________________
 509 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
 510 {
 511   //
 512   // Get First two daughters of the mother
 513   //
 514
 515   Int_t lblD1=-1;
 516   Int_t lblD2=-1;
 517   d1=0;
 518   d2=0;
 519   if (!fMCEvent) return;
 520   if (fAnaType==kAOD){
 521     const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
 522     lblD1=aodMother->GetDaughter(0);
 523     lblD2=aodMother->GetDaughter(1);
 524   } else if (fAnaType==kESD){
 525     const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
 526     lblD1=aodMother->GetFirstDaughter();
 527     lblD2=aodMother->GetLastDaughter();
 528   }
 529   d1=fMCEvent->GetTrack(lblD1);
 530   d2=fMCEvent->GetTrack(lblD2);
 531 }
 532
 533 //____________________________________________________________
 534 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair)
 535 {
 536   //
 537   // Check whether two particles have the same mother
 538   //
 539
 540   const AliVParticle * daughter1 = pair->GetFirstDaughter();
 541   const AliVParticle * daughter2 = pair->GetSecondDaughter();
 542
 543   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
 544   AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
 545   if (!mcDaughter1 || !mcDaughter2) return 0;
 546
 547   Int_t labelMother1=-1;
 548   Int_t labelMother2=-1;
 549
 550   if (mcDaughter1->IsA()==AliMCParticle::Class()){
 551     labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
 552     labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
 553   } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
 554     labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
 555     labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
 556   }
 557
 558   Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
 559
 560   return sameMother;
 561 }
 562
 563
 564
 565
 566
 567
 568
 569