PWGDQ/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 <AliAODInputHandler.h>
  29 #include <AliMCEventHandler.h>
  30 #include <AliMCEvent.h>
  31 #include <AliMCParticle.h>
  32 #include <AliAODMCParticle.h>
  33 #include <AliAODMCHeader.h>
  34 #include <AliStack.h>
  35 #include <AliESDEvent.h>
  36 #include <AliAODEvent.h>
  37 #include <AliESDtrack.h>
  38 #include <AliAODTrack.h>
  39 #include <AliLog.h>
  40
  41 #include <TClonesArray.h>
  42 #include <TParticle.h>
  43
  44 #include "AliDielectronSignalMC.h"
  45 #include "AliDielectronMC.h"
  46
  47 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
  48
  49 //____________________________________________________________
  50 AliDielectronMC* AliDielectronMC::Instance()
  51 {
  52   //
  53   // return pointer to singleton implementation
  54   //
  55   if (fgInstance) return fgInstance;
  56
  57   AnalysisType type=kUNSET;
  58   Bool_t hasMC=kFALSE;
  59   if (AliAnalysisManager::GetAnalysisManager()){
  60     if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
  61     else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
  62
  63     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  64     if(type == kESD) hasMC=mcHandler!=0x0;
  65     }
  66
  67   fgInstance=new AliDielectronMC(type);
  68
  69   fgInstance->SetHasMC(hasMC);
  70
  71   return fgInstance;
  72 }
  73
  74 //____________________________________________________________
  75 AliDielectronMC::AliDielectronMC(AnalysisType type):
  76   fMCEvent(0x0),
  77   fStack(0x0),
  78   fAnaType(type),
  79   fHasMC(kTRUE),
  80   fMcArray(0x0)
  81 {
  82   //
  83   // default constructor
  84   //
  85 }
  86
  87
  88 //____________________________________________________________
  89 AliDielectronMC::~AliDielectronMC()
  90 {
  91   //
  92   // default destructor
  93   //
  94
  95 }
  96
  97 //____________________________________________________________
  98 void AliDielectronMC::Initialize()
  99 {
 100   //
 101   // initialize MC class
 102   //
 103   if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
 104 }
 105
 106 //____________________________________________________________
 107 Int_t AliDielectronMC::GetNMCTracks()
 108 {
 109   //
 110   //  return the number of generated tracks from MC event
 111   //
 112   if(fAnaType == kESD){
 113     if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
 114     return fMCEvent->GetNumberOfTracks();}
 115   else if(fAnaType == kAOD){
 116     if(!fMcArray) { AliError("No fMcArray"); return 0; }
 117     return fMcArray->GetEntriesFast();
 118   }
 119   return 0;
 120 }
 121
 122 //____________________________________________________________
 123 Int_t AliDielectronMC::GetNMCTracksFromStack()
 124 {
 125   //
 126   //  return the number of generated tracks from stack
 127   //
 128   if (!fStack){ AliError("No fStack"); return -999; }
 129   return fStack->GetNtrack();
 130 }
 131
 132 //____________________________________________________________
 133 Int_t AliDielectronMC::GetNPrimary() const
 134 {
 135   //
 136   //  return the number of primary track from MC event
 137   //
 138   if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
 139   return fMCEvent->GetNumberOfPrimaries();
 140 }
 141
 142 //____________________________________________________________
 143 Int_t AliDielectronMC::GetNPrimaryFromStack()
 144 {
 145   //
 146   //  return the number of primary track from stack
 147   //
 148   if (!fStack){ AliError("No fStack"); return -999; }
 149   return fStack->GetNprimary();
 150 }
 151
 152 //____________________________________________________________
 153 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk) const
 154 {
 155   //
 156   // return MC track directly from MC event
 157   //
 158   if (itrk<0) return NULL;
 159   AliVParticle * track=0x0;
 160   if(fAnaType == kESD){
 161     if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 162     track = fMCEvent->GetTrack(itrk); //  tracks from MC event (ESD)
 163   } else if(fAnaType == kAOD) {
 164     if (!fMcArray){ AliError("No fMcArray"); return NULL;}
 165     if (itrk>fMcArray->GetEntriesFast()) { AliWarning(Form("track %d out of array size %d",itrk,fMcArray->GetEntriesFast())); return NULL;}
 166     track = (AliVParticle*)fMcArray->At(itrk); //  tracks from MC event (AOD)
 167   }
 168   return track;
 169 }
 170
 171 //____________________________________________________________
 172 Bool_t AliDielectronMC::ConnectMCEvent()
 173 {
 174   //
 175   // connect stack object from the mc handler
 176   //
 177
 178   fMcArray = 0x0;
 179   fMCEvent = 0x0;
 180
 181   if(fAnaType == kESD){
 182     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 183     if (!mcHandler){ /*AliError("Could not retrive MC event handler!");*/ return kFALSE; }
 184     if (!mcHandler->InitOk() ) return kFALSE;
 185     if (!mcHandler->TreeK() )  return kFALSE;
 186     if (!mcHandler->TreeTR() ) return kFALSE;
 187
 188     AliMCEvent* mcEvent = mcHandler->MCEvent();
 189     if (!mcEvent){ /*AliError("Could not retrieve MC event!");*/ return kFALSE; }
 190     fMCEvent = mcEvent;
 191
 192     if (!UpdateStack()) return kFALSE;
 193   }
 194   else if(fAnaType == kAOD)
 195   {
 196     AliAODInputHandler* aodHandler=(AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 197     if (!aodHandler) return kFALSE;
 198     AliAODEvent *aod=aodHandler->GetEvent();
 199     if (!aod) return kFALSE;
 200
 201     fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 202     if (!fMcArray){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
 203     else fHasMC=kTRUE;
 204   }
 205   return kTRUE;
 206 }
 207
 208 //____________________________________________________________
 209 Bool_t AliDielectronMC::UpdateStack()
 210 {
 211   //
 212   // update stack with new event
 213   //
 214   if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
 215   AliStack* stack = fMCEvent->Stack();
 216   if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
 217   fStack = stack;
 218   return kTRUE;
 219 }
 220
 221 //____________________________________________________________
 222 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
 223 {
 224   //
 225   // return MC track
 226   //
 227   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 228
 229   Int_t nStack = fMCEvent->GetNumberOfTracks();
 230   Int_t label = TMath::Abs(_track->GetLabel());
 231   if(label>nStack)return NULL;
 232   if(label<0) return NULL;
 233   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
 234   return mctrack;
 235 }
 236
 237
 238 //____________________________________________________________
 239 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
 240 {
 241   //
 242   // return MC track
 243   //
 244  if(!fMcArray) { AliError("No fMcArray"); return NULL;}
 245  Int_t label = _track->GetLabel();
 246  if(label < 0 || label > fMcArray->GetEntriesFast()) return NULL;
 247  AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
 248  return mctrack;
 249 }
 250
 251 //____________________________________________________________
 252 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
 253 {
 254   //
 255   // return MC track from stack
 256   //
 257   Int_t label = TMath::Abs(_track->GetLabel());
 258   if (!fStack) AliWarning("fStack is not available. Update stack first.");
 259   TParticle* mcpart = fStack->Particle(label);
 260   if (!mcpart) return NULL;
 261   return mcpart;
 262 }
 263
 264 //____________________________________________________________
 265 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
 266 {
 267   //
 268   // return MC track mother
 269   //
 270   AliMCParticle* mcpart = GetMCTrack(_track);
 271   if (!mcpart) return NULL;
 272   if(mcpart->GetMother()<0) return NULL;
 273   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
 274   if (!mcmother) return NULL;
 275   return mcmother;
 276 }
 277
 278 //______________________________________________________________
 279 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
 280 {
 281  //
 282  // return MC track mother
 283  //
 284  AliAODMCParticle* mcpart = GetMCTrack(_track);
 285  if (!mcpart) return NULL;
 286  if(mcpart->GetMother() < 0) return NULL;
 287  AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
 288  if (!mcmother) return NULL;
 289  return mcmother;
 290 }
 291 //____________________________________________________________
 292 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
 293   //
 294   // return MC track mother
 295   //
 296   if(_particle->GetMother() < 0) return NULL;
 297   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
 298   return mcmother;
 299 }
 300
 301 //____________________________________________________________
 302 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
 303   //
 304   // return MC track mother
 305   //
 306   if( _particle->GetMother() < 0) return NULL;
 307   AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
 308   return mcmother;
 309 }
 310
 311 //____________________________________________________________
 312 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
 313 {
 314   //
 315   // return MC track mother from stack
 316   //
 317   TParticle* mcpart = GetMCTrackFromStack(_track);
 318   if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL;
 319   TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
 320   if (!mcmother) return NULL;
 321   return mcmother;
 322 }
 323
 324 //____________________________________________________________
 325 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
 326 {
 327   //
 328   // return PDG code of the track from the MC truth info
 329   //
 330   AliMCParticle* mcpart = GetMCTrack(_track);
 331   if (!mcpart) return -999;
 332   return mcpart->PdgCode();
 333 }
 334
 335 //__________________________________________________________
 336 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
 337 {
 338  //
 339  // return PDG code of the track from the MC truth info
 340  //
 341  AliAODMCParticle* mcpart = GetMCTrack(_track);
 342  if (!mcpart) return -999;
 343  return mcpart->PdgCode();
 344 }
 345
 346 //____________________________________________________________
 347 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
 348 {
 349   //
 350   // return MC PDG code from stack
 351   //
 352   TParticle* mcpart = GetMCTrackFromStack(_track);
 353   if (!mcpart) return -999;
 354   return mcpart->GetPdgCode();
 355 }
 356
 357 //____________________________________________________________
 358 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
 359 {
 360   //
 361   // return PDG code of the mother track from the MC truth info
 362   //
 363   AliMCParticle* mcmother = GetMCTrackMother(_track);
 364   if (!mcmother) return -999;
 365   return mcmother->PdgCode();
 366 }
 367
 368 //________________________________________________________
 369 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
 370 {
 371   //
 372   // return PDG code of the mother track from the MC truth info
 373   //
 374   AliAODMCParticle* mcmother = GetMCTrackMother(_track);
 375   if (!mcmother) return -999;
 376   return mcmother->PdgCode();
 377 }
 378
 379 //____________________________________________________________
 380 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
 381 {
 382   //
 383   // return PDG code of the mother track from stack
 384   //
 385   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 386   if (!mcmother) return -999;
 387   return mcmother->GetPdgCode();
 388 }
 389
 390 //____________________________________________________________
 391 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
 392 {
 393   //
 394   // return process number of the track
 395   //
 396   AliMCParticle* mcpart = GetMCTrack(_track);
 397   if (!mcpart) return -999;
 398   return 0;
 399 }
 400
 401 //____________________________________________________________
 402 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
 403 {
 404   //
 405   // return process number of the track
 406   //
 407   TParticle* mcpart = GetMCTrackFromStack(_track);
 408   if (!mcpart) return -999;
 409   return mcpart->GetUniqueID();
 410 }
 411
 412 //____________________________________________________________
 413 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
 414 {
 415   //
 416   // returns the number of daughters
 417   //
 418   AliMCParticle *mcmother=GetMCTrackMother(track);
 419   if(!mcmother||!mcmother->Particle()) return -999;
 420 //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 421   return mcmother->Particle()->GetNDaughters();
 422 }
 423
 424 //_________________________________________________________
 425 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
 426 {
 427   //
 428   // returns the number of daughters
 429   //
 430   AliAODMCParticle *mcmother=GetMCTrackMother(track);
 431   if(!mcmother) return -999;
 432   return NumberOfDaughters(mcmother);
 433
 434 }
 435
 436 //____________________________________________________________
 437 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
 438 {
 439   //
 440   // returns the number of daughters
 441   //
 442   AliMCParticle *mcmother=GetMCTrackMother(particle);
 443   if(!mcmother||!mcmother->Particle()) return -999;
 444   //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 445   return mcmother->Particle()->GetNDaughters();
 446 }
 447
 448 //____________________________________________________________
 449 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
 450 {
 451   //
 452   // returns the number of daughters
 453   //
 454   AliAODMCParticle *mcmother=GetMCTrackMother(particle);
 455   if(!mcmother) return -999;
 456   return mcmother->GetNDaughters();
 457 }
 458
 459 //____________________________________________________________
 460 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
 461 {
 462   //
 463   // return process number of the mother of the track
 464   //
 465   AliMCParticle* mcmother = GetMCTrackMother(_track);
 466   if (!mcmother) return -999;
 467   return 0;
 468 }
 469
 470 //____________________________________________________________
 471 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
 472 {
 473   //
 474   // return process number of the mother of the track
 475   //
 476   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 477   if (!mcmother) return -999;
 478   return mcmother->GetUniqueID();
 479 }
 480
 481 //____________________________________________________________
 482 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
 483 {
 484   //
 485   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 486   //
 487   if (fAnaType==kESD && !fMCEvent) return kFALSE;
 488   if (fAnaType==kAOD && !fMcArray) return kFALSE;
 489   if (!particle) return kFALSE;
 490
 491   if (particle->IsA()==AliMCParticle::Class()){
 492     return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
 493   } else if (particle->IsA()==AliAODMCParticle::Class()){
 494    return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
 495   } else {
 496     AliError("Unknown particle type");
 497   }
 498   return kFALSE;
 499 }
 500
 501 //____________________________________________________________
 502 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
 503 {
 504   //
 505   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 506   // ESD case
 507   //
 508
 509   //check pdg code
 510   if (particle->PdgCode()!=pdgMother) return kFALSE;
 511   Int_t ifirst = particle->GetFirstDaughter();
 512   Int_t ilast  = particle->GetLastDaughter();
 513
 514   //check number of daughters
 515   if ((ilast-ifirst)!=1) return kFALSE;
 516   AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 517   AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
 518
 519   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 520   if (firstD->Charge()>0){
 521     if (firstD->PdgCode()!=-11) return kFALSE;
 522     if (secondD->PdgCode()!=11) return kFALSE;
 523   }else{
 524     if (firstD->PdgCode()!=11) return kFALSE;
 525     if (secondD->PdgCode()!=-11) return kFALSE;
 526   }
 527
 528   return kTRUE;
 529 }
 530
 531 //____________________________________________________________
 532 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
 533 {
 534   //
 535   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 536   // AOD case
 537   //
 538
 539   if (particle->GetPdgCode()!=pdgMother) return kFALSE;
 540   if (particle->GetNDaughters()!=2) return kFALSE;
 541
 542   Int_t ifirst = particle->GetDaughter(0);
 543   Int_t ilast  = particle->GetDaughter(1);
 544
 545   //check number of daughters
 546   if ((ilast-ifirst)!=1) return kFALSE;
 547
 548   AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 549   AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
 550
 551   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 552
 553   if (firstD->Charge()>0){
 554     if (firstD->GetPdgCode()!=-11) return kFALSE;
 555     if (secondD->GetPdgCode()!=11) return kFALSE;
 556   }else{
 557     if (firstD->GetPdgCode()!=11) return kFALSE;
 558     if (secondD->GetPdgCode()!=-11) return kFALSE;
 559   }
 560   return kTRUE;
 561 }
 562
 563 //____________________________________________________________
 564 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 565 {
 566   //
 567   // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
 568   //
 569   if (fAnaType==kESD){
 570   if (!fMCEvent) return -1;
 571   return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
 572   }
 573   else if (fAnaType==kAOD)
 574   {
 575   if (!fMcArray) return -1;
 576   return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
 577   }
 578
 579   return -1;
 580 }
 581
 582 //____________________________________________________________
 583 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 584 {
 585   //
 586   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 587   //    have the same mother and the mother had pdg code pdgMother
 588   // ESD case
 589   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 590   //
 591
 592   AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
 593   AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
 594
 595   if (!mcPart1||!mcPart2) return -1;
 596
 597   Int_t lblMother1=mcPart1->GetMother();
 598   Int_t lblMother2=mcPart2->GetMother();
 599
 600   AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 601   if (!mcMother1) return -1;
 602   if (lblMother1!=lblMother2) return -1;
 603   if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
 604   if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
 605   if (mcMother1->PdgCode()!=pdgMother) return -1;
 606
 607   return lblMother1;
 608 }
 609
 610 //____________________________________________________________
 611 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 612 {
 613   //
 614   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 615   //    have the same mother and the mother had pdg code pdgMother
 616   // AOD case
 617   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 618   //
 619   AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
 620   AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
 621
 622   if (!mcPart1||!mcPart2) return -1;
 623
 624   Int_t lblMother1=mcPart1->GetMother();
 625   Int_t lblMother2=mcPart2->GetMother();
 626
 627   AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 628
 629   if (!mcMother1) return -1;
 630   if (lblMother1!=lblMother2) return -1;
 631   if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
 632   if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
 633   if (mcMother1->GetPdgCode()!=pdgMother) return -1;
 634
 635   return lblMother1;
 636 }
 637
 638 //____________________________________________________________
 639 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
 640 {
 641   //
 642   // Get First two daughters of the mother
 643   //
 644   Int_t lblD1=-1;
 645   Int_t lblD2=-1;
 646   d1=0;
 647   d2=0;
 648   if (fAnaType==kAOD){
 649     if(!fMcArray) return;
 650     const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
 651     lblD1=aodMother->GetDaughter(0);
 652     lblD2=aodMother->GetDaughter(1);
 653     d1 = (AliVParticle*)fMcArray->At(lblD1);
 654     d2 = (AliVParticle*)fMcArray->At(lblD2);
 655    } else if (fAnaType==kESD){
 656     if (!fMCEvent) return;
 657     const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
 658     lblD1=aodMother->GetFirstDaughter();
 659     lblD2=aodMother->GetLastDaughter();
 660     d1=fMCEvent->GetTrack(lblD1);
 661     d2=fMCEvent->GetTrack(lblD2);
 662    }
 663 }
 664
 665
 666 //________________________________________________________________________________
 667 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
 668   //
 669   //  Get the label of the mother for particle with label daughterLabel
 670   //
 671   if(daughterLabel<0) return -1;
 672   if (fAnaType==kAOD) {
 673     if(!fMcArray) return -1;
 674     if(GetMCTrackFromMCEvent(daughterLabel))
 675       return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
 676   } else if(fAnaType==kESD) {
 677     if (!fMCEvent) return -1;
 678     if(GetMCTrackFromMCEvent(daughterLabel))
 679       return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
 680   }
 681   return -1;
 682 }
 683
 684
 685 //________________________________________________________________________________
 686 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
 687   //
 688   //  Get particle code using the label from stack
 689   //
 690   if(label<0) return 0;
 691   if(fAnaType==kAOD) {
 692     if(!fMcArray) return 0;
 693     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
 694   } else if(fAnaType==kESD) {
 695     if (!fMCEvent) return 0;
 696     return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
 697   }
 698   return 0;
 699 }
 700
 701
 702 //________________________________________________________________________________
 703 Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t pdgExclusion, Bool_t checkBothCharges) const {
 704   //
 705   //  Test the PDG codes of particles with the required ones
 706   //
 707   Bool_t result = kTRUE;
 708   Int_t absRequiredPDG = TMath::Abs(requiredPDG);
 709
 710   switch(absRequiredPDG) {
 711   case 0:
 712     result = kTRUE;    // PDG not required (any code will do fine)
 713     break;
 714   case 100:     // light flavoured mesons
 715     if(checkBothCharges)
 716       result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=199;
 717     else {
 718       if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=199;
 719       if(requiredPDG<0) result = particlePDG>=-199 && particlePDG<=-100;
 720     }
 721     break;
 722   case 1000:     // light flavoured baryons
 723     if(checkBothCharges)
 724       result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=1999;
 725     else {
 726       if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=1999;
 727       if(requiredPDG<0) result = particlePDG>=-1999 && particlePDG<=-1000;
 728     }
 729     break;
 730   case 200:     // light flavoured mesons
 731     if(checkBothCharges)
 732       result = TMath::Abs(particlePDG)>=200 && TMath::Abs(particlePDG)<=299;
 733     else {
 734       if(requiredPDG>0)result = particlePDG>=200 && particlePDG<=299;
 735       if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-200;
 736     }
 737     break;
 738   case 2000:     // light flavoured baryons
 739     if(checkBothCharges)
 740       result = TMath::Abs(particlePDG)>=2000 && TMath::Abs(particlePDG)<=2999;
 741     else {
 742       if(requiredPDG>0) result = particlePDG>=2000 && particlePDG<=2999;
 743       if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-2000;
 744     }
 745     break;
 746   case 300:     // all strange mesons
 747     if(checkBothCharges)
 748       result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
 749     else {
 750       if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
 751       if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
 752     }
 753     break;
 754   case 3000:     // all strange baryons
 755     if(checkBothCharges)
 756       result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
 757     else {
 758       if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
 759       if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
 760     }
 761     break;
 762   case 400:     // all charmed mesons
 763     if(checkBothCharges)
 764       result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
 765     else {
 766       if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
 767       if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
 768     }
 769     break;
 770   case 401:     // open charm mesons
 771     if(checkBothCharges)
 772       result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439;
 773     else {
 774       if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=439;
 775       if(requiredPDG<0) result = particlePDG>=-439 && particlePDG<=-400;
 776     }
 777     break;
 778   case 402:     // open charm mesons and baryons together
 779     if(checkBothCharges)
 780       result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
 781                (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399);
 782     else {
 783       if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
 784                                  (particlePDG>=4000 && particlePDG<=4399);
 785       if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
 786                                  (particlePDG>=-4399 && particlePDG<=-4000);
 787     }
 788     break;
 789   case 403:     // all charm hadrons
 790     if(checkBothCharges)
 791       result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499) ||
 792                (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999);
 793     else {
 794       if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=499) ||
 795                                  (particlePDG>=4000 && particlePDG<=4999);
 796       if(requiredPDG<0) result = (particlePDG>=-499 && particlePDG<=-400) ||
 797                                  (particlePDG>=-4999 && particlePDG<=-4000);
 798     }
 799     break;
 800   case 4000:     // all charmed baryons
 801     if(checkBothCharges)
 802       result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
 803     else {
 804       if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
 805       if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
 806     }
 807     break;
 808   case 4001:     // open charm baryons
 809     if(checkBothCharges)
 810       result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399;
 811     else {
 812       if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4399;
 813       if(requiredPDG<0) result = particlePDG>=-4399 && particlePDG<=-4000;
 814     }
 815     break;
 816   case 500:      // all beauty mesons
 817     if(checkBothCharges)
 818       result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
 819     else {
 820       if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
 821       if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
 822     }
 823     break;
 824   case 501:      // open beauty mesons
 825     if(checkBothCharges)
 826       result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549;
 827     else {
 828       if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=549;
 829       if(requiredPDG<0) result = particlePDG>=-549 && particlePDG<=-500;
 830     }
 831     break;
 832   case 502:      // open beauty mesons and baryons
 833     if(checkBothCharges)
 834       result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
 835                (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
 836     else {
 837       if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=549) ||
 838                                  (particlePDG>=5000 && particlePDG<=5499);
 839       if(requiredPDG<0) result = (particlePDG>=-549 && particlePDG<=-500) ||
 840                                  (particlePDG>=-5499 && particlePDG<=-5000);
 841     }
 842     break;
 843   case 503:      // all beauty hadrons
 844     if(checkBothCharges)
 845       result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599) ||
 846                (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999);
 847     else {
 848       if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=599) ||
 849                                  (particlePDG>=5000 && particlePDG<=5999);
 850       if(requiredPDG<0) result = (particlePDG>=-599 && particlePDG<=-500) ||
 851                                  (particlePDG>=-5999 && particlePDG<=-5000);
 852     }
 853     break;
 854   case 5000:      // all beauty baryons
 855     if(checkBothCharges)
 856       result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
 857     else {
 858       if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
 859       if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
 860     }
 861     break;
 862   case 5001:      // open beauty baryons
 863     if(checkBothCharges)
 864       result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499;
 865     else {
 866       if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5499;
 867       if(requiredPDG<0) result = particlePDG>=-5499 && particlePDG<=-5000;
 868     }
 869     break;
 870   default:          // all specific cases
 871     if(checkBothCharges)
 872       result = (absRequiredPDG==TMath::Abs(particlePDG));
 873     else
 874       result = (requiredPDG==particlePDG);
 875   }
 876
 877   if(absRequiredPDG!=0 && pdgExclusion) result = !result;
 878   return result;
 879 }
 880
 881
 882 //________________________________________________________________________________
 883 Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
 884   //
 885   // Check if the particle with label "label" is a physical primary according to the
 886   // definition in AliStack::IsPhysicalPrimary(Int_t label)
 887   // Convention for being physical primary:
 888   // 1.) particles produced in the collision
 889   // 2.) stable particles with respect to strong and electromagnetic interactions
 890   // 3.) excludes initial state particles
 891   // 4.) includes products of directly produced Sigma0 hyperon decay
 892   // 5.) includes products of directly produced pi0 decays
 893   // 6.) includes products of directly produced beauty hadron decays
 894   //
 895   if(label<0) return kFALSE;
 896   if(fAnaType==kAOD) {
 897     if(!fMcArray) return kFALSE;
 898     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsPhysicalPrimary();
 899   } else if(fAnaType==kESD) {
 900     if (!fMCEvent) return kFALSE;
 901     return fStack->IsPhysicalPrimary(label);
 902   }
 903   return kFALSE;
 904 }
 905
 906
 907 //________________________________________________________________________________
 908 Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
 909   //
 910   //  Check the source for the particle
 911   //
 912
 913   switch (source) {
 914     case AliDielectronSignalMC::kDontCare :
 915       return kTRUE;
 916     break;
 917     case AliDielectronSignalMC::kPrimary :
 918       // true if label is in the list of particles from physics generator
 919       // NOTE: This includes all physics event history (initial state particles,
 920       //       exchange bosons, quarks, di-quarks, strings, un-stable particles, final state particles)
 921       //       Only the final state particles make it to the detector!!
 922       return (label>=0 && label<=GetNPrimary());
 923     break;
 924     case AliDielectronSignalMC::kFinalState :
 925       // primary particles created in the collision which reach the detectors
 926       // These would be:
 927       // 1.) particles produced in the collision
 928       // 2.) stable particles with respect to strong and electromagnetic interactions
 929       // 3.) excludes initial state particles
 930       // 4.) includes products of directly produced Sigma0 hyperon decay
 931       // 5.) includes products of directly produced pi0 decays
 932       // 6.) includes products of directly produced beauty hadron decays
 933       return IsPhysicalPrimary(label);
 934     break;
 935     case AliDielectronSignalMC::kDirect :
 936       // Primary particles which do not have any mother
 937       // This is the case for:
 938       // 1.) Initial state particles (the 2 protons in Pythia pp collisions)
 939       // 2.) In some codes, with sudden freeze-out, all particles generated from the fireball are direct.
 940       //     There is no history for these particles.
 941       // 3.) Certain particles added via MC generator cocktails (e.g. J/psi added to pythia MB events)
 942       return (label>=0 && GetMothersLabel(label)<0);
 943       break;
 944     case AliDielectronSignalMC::kNoCocktail :
 945       // Particles from the HIJING event and NOT from the AliGenCocktail
 946       return (label>=0 && GetMothersLabel(label)>=0);
 947       break;
 948     case AliDielectronSignalMC::kSecondary :
 949       // particles which are created by the interaction of final state primaries with the detector
 950       // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
 951       return (label>=GetNPrimary() && !IsPhysicalPrimary(label));
 952     break;
 953     default :
 954       return kFALSE;
 955   }
 956   return kFALSE;
 957 }
 958
 959 //________________________________________________________________________________
 960 Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
 961 {
 962   //
 963   // Check if the particle has a three body decay, one being a photon
 964   //
 965   if(label<0) return kFALSE;
 966
 967
 968   if(fAnaType==kAOD) {
 969     if(!fMcArray) return kFALSE;
 970     AliAODMCParticle *mother=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label));
 971     if (!mother) return kFALSE;
 972     const Int_t nd=mother->GetNDaughters();
 973     if (nd==2) return kFALSE;
 974     for (Int_t i=2; i<nd; ++i)
 975       if (GetMCTrackFromMCEvent(mother->GetDaughter(0)+i)->PdgCode()!=22) return kFALSE; //last daughter is photon
 976   } else if(fAnaType==kESD) {
 977     if (!fMCEvent) return kFALSE;
 978     AliMCParticle *mother=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label));
 979     if (!mother) return kFALSE;
 980     const Int_t nd=(mother->GetLastDaughter()-mother->GetFirstDaughter()+1);
 981     if (nd==2) return kFALSE;
 982     for (Int_t i=2; i<nd; ++i)
 983       if (GetMCTrackFromMCEvent(mother->GetFirstDaughter()+i)->PdgCode()!=22) return kFALSE; //last daughters are photons
 984   }
 985   return kTRUE;
 986 }
 987
 988 //________________________________________________________________________________
 989 Bool_t AliDielectronMC::CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const
 990 {
 991   //
 992   // Check for the decision of the radiative type request
 993   //
 994
 995   if (!signalMC) return kFALSE;
 996
 997   if (signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kAll) return kTRUE;
 998
 999   Bool_t isRadiative=CheckIsRadiative(mLabel);
1000   if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsRadiative) && !isRadiative) return kFALSE;
1001   if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsNotRadiative) && isRadiative) return kFALSE;
1002
1003   return kTRUE;
1004 }
1005
1006 //________________________________________________________________________________
1007 Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
1008   //
1009   // Check if the particle corresponds to the MC truth in signalMC in the branch specified
1010   //
1011
1012   // NOTE:  Some particles have the sign of the label flipped. It is related to the quality of matching
1013   //        between the ESD and the MC track. The negative labels indicate a poor matching quality
1014   //if(label<0) return kFALSE;
1015   if(label<0) label *= -1;
1016
1017   AliVParticle* part = GetMCTrackFromMCEvent(label);
1018   if (!part) {
1019     AliError(Form("Could not find MC particle with label %d",label));
1020     return kFALSE;
1021   }
1022
1023   // check the leg
1024   if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
1025   if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
1026
1027   // check the mother
1028   AliVParticle* mcMother=0x0;
1029   Int_t mLabel = -1;
1030   if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1031     if(part) {
1032       mLabel = GetMothersLabel(label);
1033       mcMother = GetMCTrackFromMCEvent(mLabel);
1034     }
1035     if(!mcMother && !signalMC->GetMotherPDGexclude(branch)) return kFALSE;
1036
1037     if(!ComparePDG((mcMother ? mcMother->PdgCode() : 0),signalMC->GetMotherPDG(branch),signalMC->GetMotherPDGexclude(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
1038     if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
1039
1040     //check for radiative deday
1041     if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
1042   }
1043
1044   // check the grandmother
1045   AliVParticle* mcGrandMother=0x0;
1046   if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1047     Int_t gmLabel = -1;
1048     if(mcMother) {
1049       gmLabel = GetMothersLabel(mLabel);
1050       mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
1051     }
1052     if(!mcGrandMother && !signalMC->GetGrandMotherPDGexclude(branch)) return kFALSE;
1053
1054     if(!ComparePDG((mcGrandMother ? mcGrandMother->PdgCode() : 0),signalMC->GetGrandMotherPDG(branch),signalMC->GetGrandMotherPDGexclude(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
1055     if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
1056   }
1057
1058   return kTRUE;
1059 }
1060
1061
1062
1063 //________________________________________________________________________________
1064 Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
1065   //
1066   // Check if the pair corresponds to the MC truth in signalMC
1067   //
1068
1069   // legs (daughters)
1070   const AliVParticle * mcD1 = pair->GetFirstDaughter();
1071   const AliVParticle * mcD2 = pair->GetSecondDaughter();
1072   Int_t labelD1 = (mcD1 ? mcD1->GetLabel() : -1);
1073   Int_t labelD2 = (mcD2 ? mcD2->GetLabel() : -1);
1074   if(labelD1<0) labelD1 *= -1;
1075   if(labelD2<0) labelD2 *= -1;
1076   Int_t d1Pdg = 0;
1077   Int_t d2Pdg = 0;
1078   d1Pdg=GetPdgFromLabel(labelD1);
1079   d2Pdg=GetPdgFromLabel(labelD2);
1080
1081   // mothers
1082   AliVParticle* mcM1=0x0;
1083   AliVParticle* mcM2=0x0;
1084
1085   // grand-mothers
1086   AliVParticle* mcG1 = 0x0;
1087   AliVParticle* mcG2 = 0x0;
1088
1089   // make direct(1-1 and 2-2) and cross(1-2 and 2-1) comparisons for the whole branch
1090   Bool_t directTerm = kTRUE;
1091   // daughters
1092   directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1093                && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
1094
1095   directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1096                && CheckParticleSource(labelD2, signalMC->GetLegSource(2));
1097
1098   // mothers
1099   Int_t labelM1 = -1;
1100   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1101     labelM1 = GetMothersLabel(labelD1);
1102     if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1103     directTerm = directTerm && (mcM1 || signalMC->GetMotherPDGexclude(1))
1104                  && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1105                  && CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
1106                  && CheckRadiativeDecision(labelM1,signalMC);
1107   }
1108
1109   Int_t labelM2 = -1;
1110   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1111     labelM2 = GetMothersLabel(labelD2);
1112     if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1113     directTerm = directTerm && (mcM2 || signalMC->GetMotherPDGexclude(2))
1114                  && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1115                  && CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
1116                  && CheckRadiativeDecision(labelM2,signalMC);
1117   }
1118
1119   // grand-mothers
1120   Int_t labelG1 = -1;
1121   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1122     labelG1 = GetMothersLabel(labelM1);
1123     if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1124     directTerm = directTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(1))
1125                  && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1126                  && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
1127   }
1128
1129   Int_t labelG2 = -1;
1130   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1131     labelG2 = GetMothersLabel(labelM2);
1132     if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1133     directTerm = directTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(2))
1134                  && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1135                  && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
1136   }
1137
1138   // Cross term
1139   Bool_t crossTerm = kTRUE;
1140   // daughters
1141   crossTerm = crossTerm && mcD2
1142               && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1143               && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
1144
1145   crossTerm = crossTerm && mcD1
1146               && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1147               && CheckParticleSource(labelD1, signalMC->GetLegSource(2));
1148
1149   // mothers
1150   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1151     if(!mcM2 && labelD2>-1) {
1152       labelM2 = GetMothersLabel(labelD2);
1153       if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1154     }
1155     crossTerm = crossTerm && (mcM2 || signalMC->GetMotherPDGexclude(1))
1156                 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1157                 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
1158                 && CheckRadiativeDecision(labelM2,signalMC);
1159   }
1160
1161   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1162     if(!mcM1 && labelD1>-1) {
1163       labelM1 = GetMothersLabel(labelD1);
1164       if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1165     }
1166     crossTerm = crossTerm && (mcM1 || signalMC->GetMotherPDGexclude(2))
1167                 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1168                 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
1169                 && CheckRadiativeDecision(labelM1,signalMC);
1170   }
1171
1172   // grand-mothers
1173   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1174     if(!mcG2 && mcM2) {
1175       labelG2 = GetMothersLabel(labelM2);
1176       if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1177     }
1178     crossTerm = crossTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(1))
1179                 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1180                 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
1181   }
1182
1183   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1184     if(!mcG1 && mcM1) {
1185       labelG1 = GetMothersLabel(labelM1);
1186       if(labelG2>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1187     }
1188     crossTerm = crossTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(2))
1189                 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1190                 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
1191   }
1192
1193   Bool_t motherRelation = kTRUE;
1194   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
1195     motherRelation = motherRelation && HaveSameMother(pair);
1196   }
1197   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
1198     motherRelation = motherRelation && !HaveSameMother(pair);
1199   }
1200
1201   return ((directTerm || crossTerm) && motherRelation);
1202 }
1203
1204
1205
1206 //____________________________________________________________
1207 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
1208 {
1209   //
1210   // Check whether two particles have the same mother
1211   //
1212
1213   const AliVParticle * daughter1 = pair->GetFirstDaughter();
1214   const AliVParticle * daughter2 = pair->GetSecondDaughter();
1215   if (!daughter1 || !daughter2) return 0;
1216
1217   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1218   AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
1219   if (!mcDaughter1 || !mcDaughter2) return 0;
1220
1221   Int_t labelMother1=-1;
1222   Int_t labelMother2=-1;
1223
1224   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1225     labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1226     labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
1227   } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1228     labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1229     labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
1230   }
1231
1232   Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
1233
1234   return sameMother;
1235 }
1236
1237 //________________________________________________________________
1238 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
1239 {
1240  // return: "0" for primary jpsi
1241  //         "1" for secondary jpsi (from beauty)
1242  //         "2" for background
1243  if(!HaveSameMother(pair)) return 2;
1244  AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughter())->GetLabel());
1245  Int_t labelMother=-1;
1246
1247   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1248      labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1249      } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1250      labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1251      }
1252
1253  AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
1254  if(!IsMCMotherToEE(mcMother,443)) return 2;
1255  return IsJpsiPrimary(mcMother);
1256 }
1257
1258 //______________________________________________________________
1259 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
1260 {
1261   // return: "0" for primary jpsi
1262   //         "1" for secondary jpsi (come from B decay)
1263  Int_t labelMoth=-1;
1264  Int_t pdgCode;
1265
1266  if (particle->IsA()==AliMCParticle::Class()){
1267      labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1268      while(labelMoth>0){
1269        particle = GetMCTrackFromMCEvent(labelMoth);
1270        pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
1271        if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1272        labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1273        }
1274     }
1275  else if (particle->IsA()==AliAODMCParticle::Class()){
1276      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1277      while(labelMoth>0){
1278      particle = GetMCTrackFromMCEvent(labelMoth);
1279      pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
1280      if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1281      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1282      }
1283   }
1284   return 0;
1285 }
1286
1287
1288 Bool_t AliDielectronMC::GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ){
1289
1290      if(fAnaType == kESD){
1291      const AliVVertex* mcVtx =  fMCEvent->GetPrimaryVertex();
1292      if(!mcVtx) return kFALSE;
1293      primVtxX = mcVtx->GetX();
1294      primVtxY = mcVtx->GetY();
1295      primVtxZ = mcVtx->GetZ();
1296      }else if(fAnaType == kAOD){
1297      AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
1298      if(!aod) return kFALSE;
1299      AliAODMCHeader *mcHead = dynamic_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1300      if(!mcHead) return kFALSE;
1301      primVtxX = mcHead->GetVtxX();
1302      primVtxY = mcHead->GetVtxY();
1303      primVtxZ = mcHead->GetVtxZ();
1304      }
1305
1306 return kTRUE;
1307 }