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