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