1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //#####################################################
18 //# Class AliDielectronMC #
19 //# Cut Class for Jpsi->e+e- analysis #
21 //# by WooJin J. Park, GSI / W.J.Park@gsi.de #
23 //#####################################################
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>
34 #include <AliESDEvent.h>
35 #include <AliAODEvent.h>
36 #include <AliESDtrack.h>
37 #include <AliAODTrack.h>
40 #include <TClonesArray.h>
42 #include "AliDielectronSignalMC.h"
43 #include "AliDielectronMC.h"
45 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
47 //____________________________________________________________
48 AliDielectronMC* AliDielectronMC::Instance()
51 // return pointer to singleton implementation
53 if (fgInstance) return fgInstance;
55 AnalysisType type=kUNSET;
57 if (AliAnalysisManager::GetAnalysisManager()){
58 if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
59 else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
61 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
62 if(type == kESD) hasMC=mcHandler!=0x0;
65 fgInstance=new AliDielectronMC(type);
67 fgInstance->SetHasMC(hasMC);
72 //____________________________________________________________
73 AliDielectronMC::AliDielectronMC(AnalysisType type):
81 // default constructor
86 //____________________________________________________________
87 AliDielectronMC::~AliDielectronMC()
95 //____________________________________________________________
96 void AliDielectronMC::Initialize()
99 // initialize MC class
101 if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
104 //____________________________________________________________
105 Int_t AliDielectronMC::GetNMCTracks()
108 // return the number of generated tracks from MC event
110 if(fAnaType == kESD){
111 if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
112 return fMCEvent->GetNumberOfTracks();}
113 else if(fAnaType == kAOD){
114 if(!fMcArray) { AliError("No fMcArray"); return 0; }
115 return fMcArray->GetEntriesFast();
120 //____________________________________________________________
121 Int_t AliDielectronMC::GetNMCTracksFromStack()
124 // return the number of generated tracks from stack
126 if (!fStack){ AliError("No fStack"); return -999; }
127 return fStack->GetNtrack();
130 //____________________________________________________________
131 Int_t AliDielectronMC::GetNPrimary() const
134 // return the number of primary track from MC event
136 if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
137 return fMCEvent->GetNumberOfPrimaries();
140 //____________________________________________________________
141 Int_t AliDielectronMC::GetNPrimaryFromStack()
144 // return the number of primary track from stack
146 if (!fStack){ AliError("No fStack"); return -999; }
147 return fStack->GetNprimary();
150 //____________________________________________________________
151 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk) const
154 // return MC track directly from MC event
156 if (itrk<0) return NULL;
157 AliVParticle * track=0x0;
158 if(fAnaType == kESD){
159 if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
160 track = fMCEvent->GetTrack(itrk); // tracks from MC event (ESD)
161 } else if(fAnaType == kAOD) {
162 if (!fMcArray){ AliError("No fMCEvent"); return NULL;}
163 track = (AliVParticle*)fMcArray->At(itrk); // tracks from MC event (AOD)
168 //____________________________________________________________
169 Bool_t AliDielectronMC::ConnectMCEvent()
172 // connect stack object from the mc handler
174 if(fAnaType == kESD){
176 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
177 if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
178 if (!mcHandler->InitOk() ) return kFALSE;
179 if (!mcHandler->TreeK() ) return kFALSE;
180 if (!mcHandler->TreeTR() ) return kFALSE;
182 AliMCEvent* mcEvent = mcHandler->MCEvent();
183 if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
186 if (!UpdateStack()) return kFALSE;
188 else if(fAnaType == kAOD)
191 AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
192 fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
193 if(!fMcArray) return kFALSE;
198 //____________________________________________________________
199 Bool_t AliDielectronMC::UpdateStack()
202 // update stack with new event
204 if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
205 AliStack* stack = fMCEvent->Stack();
206 if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
211 //____________________________________________________________
212 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
217 if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
219 Int_t nStack = fMCEvent->GetNumberOfTracks();
220 Int_t label = TMath::Abs(_track->GetLabel());
221 if(label>nStack)return NULL;
223 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
228 //____________________________________________________________
229 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
234 if(!fMcArray) { AliError("No fMCArray"); return NULL;}
235 Int_t label = _track->GetLabel();
236 if(label < 0) return NULL;
237 AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
241 //____________________________________________________________
242 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
245 // return MC track from stack
247 Int_t label = TMath::Abs(_track->GetLabel());
248 if (!fStack) AliWarning("fStack is not available. Update stack first.");
249 TParticle* mcpart = fStack->Particle(label);
250 if (!mcpart) return NULL;
254 //____________________________________________________________
255 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
258 // return MC track mother
260 AliMCParticle* mcpart = GetMCTrack(_track);
261 if (!mcpart) return NULL;
262 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
263 if (!mcmother) return NULL;
267 //______________________________________________________________
268 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
271 // return MC track mother
273 AliAODMCParticle* mcpart = GetMCTrack(_track);
274 if (!mcpart) return NULL;
275 if(mcpart->GetMother() < 0) return NULL;
276 AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
277 if (!mcmother) return NULL;
280 //____________________________________________________________
281 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
283 // return MC track mother
285 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
289 //____________________________________________________________
290 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
292 // return MC track mother
294 if( _particle->GetMother() < 0) return NULL;
295 AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
299 //____________________________________________________________
300 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
303 // return MC track mother from stack
305 TParticle* mcpart = GetMCTrackFromStack(_track);
306 if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL;
307 TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
308 if (!mcmother) return NULL;
312 //____________________________________________________________
313 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
316 // return PDG code of the track from the MC truth info
318 AliMCParticle* mcpart = GetMCTrack(_track);
319 if (!mcpart) return -999;
320 return mcpart->PdgCode();
323 //__________________________________________________________
324 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
327 // return PDG code of the track from the MC truth info
329 AliAODMCParticle* mcpart = GetMCTrack(_track);
330 if (!mcpart) return -999;
331 return mcpart->PdgCode();
334 //____________________________________________________________
335 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
338 // return MC PDG code from stack
340 TParticle* mcpart = GetMCTrackFromStack(_track);
341 if (!mcpart) return -999;
342 return mcpart->GetPdgCode();
345 //____________________________________________________________
346 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
349 // return PDG code of the mother track from the MC truth info
351 AliMCParticle* mcmother = GetMCTrackMother(_track);
352 if (!mcmother) return -999;
353 return mcmother->PdgCode();
356 //________________________________________________________
357 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
360 // return PDG code of the mother track from the MC truth info
362 AliAODMCParticle* mcmother = GetMCTrackMother(_track);
363 if (!mcmother) return -999;
364 return mcmother->PdgCode();
367 //____________________________________________________________
368 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
371 // return PDG code of the mother track from stack
373 TParticle* mcmother = GetMCTrackMotherFromStack(_track);
374 if (!mcmother) return -999;
375 return mcmother->GetPdgCode();
378 //____________________________________________________________
379 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
382 // return process number of the track
384 AliMCParticle* mcpart = GetMCTrack(_track);
385 if (!mcpart) return -999;
389 //____________________________________________________________
390 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
393 // return process number of the track
395 TParticle* mcpart = GetMCTrackFromStack(_track);
396 if (!mcpart) return -999;
397 return mcpart->GetUniqueID();
400 //____________________________________________________________
401 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
404 // returns the number of daughters
406 AliMCParticle *mcmother=GetMCTrackMother(track);
407 if(!mcmother||!mcmother->Particle()) return -999;
408 // return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
409 return mcmother->Particle()->GetNDaughters();
412 //_________________________________________________________
413 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
416 // returns the number of daughters
418 AliAODMCParticle *mcmother=GetMCTrackMother(track);
419 if(!mcmother) return -999;
420 return NumberOfDaughters(mcmother);
424 //____________________________________________________________
425 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
428 // returns the number of daughters
430 AliMCParticle *mcmother=GetMCTrackMother(particle);
431 if(!mcmother||!mcmother->Particle()) return -999;
432 //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
433 return mcmother->Particle()->GetNDaughters();
436 //____________________________________________________________
437 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
440 // returns the number of daughters
442 AliAODMCParticle *mcmother=GetMCTrackMother(particle);
443 if(!mcmother) return -999;
444 return mcmother->GetNDaughters();
447 //____________________________________________________________
448 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
451 // return process number of the mother of the track
453 AliMCParticle* mcmother = GetMCTrackMother(_track);
454 if (!mcmother) return -999;
458 //____________________________________________________________
459 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
462 // return process number of the mother of the track
464 TParticle* mcmother = GetMCTrackMotherFromStack(_track);
465 if (!mcmother) return -999;
466 return mcmother->GetUniqueID();
469 //____________________________________________________________
470 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
473 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
475 if (fAnaType==kESD && !fMCEvent) return kFALSE;
476 if (fAnaType==kAOD && !fMcArray) return kFALSE;
477 if (!particle) return kFALSE;
479 if (particle->IsA()==AliMCParticle::Class()){
480 return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
481 } else if (particle->IsA()==AliAODMCParticle::Class()){
482 return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
484 AliError("Unknown particle type");
489 //____________________________________________________________
490 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
493 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
498 if (particle->PdgCode()!=pdgMother) return kFALSE;
499 Int_t ifirst = particle->GetFirstDaughter();
500 Int_t ilast = particle->GetLastDaughter();
502 //check number of daughters
503 if ((ilast-ifirst)!=1) return kFALSE;
504 AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
505 AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
507 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
508 if (firstD->Charge()>0){
509 if (firstD->PdgCode()!=-11) return kFALSE;
510 if (secondD->PdgCode()!=11) return kFALSE;
512 if (firstD->PdgCode()!=11) return kFALSE;
513 if (secondD->PdgCode()!=-11) return kFALSE;
519 //____________________________________________________________
520 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
523 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
527 if (particle->GetPdgCode()!=pdgMother) return kFALSE;
528 if (particle->GetNDaughters()!=2) return kFALSE;
530 Int_t ifirst = particle->GetDaughter(0);
531 Int_t ilast = particle->GetDaughter(1);
533 //check number of daughters
534 if ((ilast-ifirst)!=1) return kFALSE;
536 AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
537 AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
539 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
541 if (firstD->Charge()>0){
542 if (firstD->GetPdgCode()!=-11) return kFALSE;
543 if (secondD->GetPdgCode()!=11) return kFALSE;
545 if (firstD->GetPdgCode()!=11) return kFALSE;
546 if (secondD->GetPdgCode()!=-11) return kFALSE;
551 //____________________________________________________________
552 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
555 // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
558 if (!fMCEvent) return -1;
559 return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
561 else if (fAnaType==kAOD)
563 if (!fMcArray) return -1;
564 return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
570 //____________________________________________________________
571 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
574 // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
575 // have the same mother and the mother had pdg code pdgMother
577 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
580 AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
581 AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
583 if (!mcPart1||!mcPart2) return -1;
585 Int_t lblMother1=mcPart1->GetMother();
586 Int_t lblMother2=mcPart2->GetMother();
588 AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
589 if (!mcMother1) return -1;
590 if (lblMother1!=lblMother2) return -1;
591 if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
592 if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
593 if (mcMother1->PdgCode()!=pdgMother) return -1;
598 //____________________________________________________________
599 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
602 // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
603 // have the same mother and the mother had pdg code pdgMother
605 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
607 AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
608 AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
610 if (!mcPart1||!mcPart2) return -1;
612 Int_t lblMother1=mcPart1->GetMother();
613 Int_t lblMother2=mcPart2->GetMother();
615 AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
617 if (!mcMother1) return -1;
618 if (lblMother1!=lblMother2) return -1;
619 if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
620 if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
621 if (mcMother1->GetPdgCode()!=pdgMother) return -1;
626 //____________________________________________________________
627 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
630 // Get First two daughters of the mother
637 if(!fMcArray) return;
638 const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
639 lblD1=aodMother->GetDaughter(0);
640 lblD2=aodMother->GetDaughter(1);
641 d1 = (AliVParticle*)fMcArray->At(lblD1);
642 d2 = (AliVParticle*)fMcArray->At(lblD2);
643 } else if (fAnaType==kESD){
644 if (!fMCEvent) return;
645 const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
646 lblD1=aodMother->GetFirstDaughter();
647 lblD2=aodMother->GetLastDaughter();
648 d1=fMCEvent->GetTrack(lblD1);
649 d2=fMCEvent->GetTrack(lblD2);
654 //________________________________________________________________________________
655 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
657 // Get the label of the mother for particle with label daughterLabel
659 if(daughterLabel<0) return -1;
660 if (fAnaType==kAOD) {
661 if(!fMcArray) return -1;
662 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
663 } else if(fAnaType==kESD) {
664 if (!fMCEvent) return -1;
665 return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
671 //________________________________________________________________________________
672 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
674 // Get particle code using the label from stack
676 if(label<0) return 0;
678 if(!fMcArray) return 0;
679 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
680 } else if(fAnaType==kESD) {
681 if (!fMCEvent) return 0;
682 return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
688 //________________________________________________________________________________
689 Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t checkBothCharges) const {
691 // Test the PDG codes of particles with the required ones
693 Bool_t result = kTRUE;
694 Int_t absRequiredPDG = TMath::Abs(requiredPDG);
695 switch(absRequiredPDG) {
697 result = kTRUE; // PDG not required (any code will do fine)
699 case 100: // all light flavoured mesons
701 result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=299;
703 if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=299;
704 if(requiredPDG<0) result = particlePDG>=-299 && particlePDG<=-100;
707 case 1000: // all light flavoured baryons
709 result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=2999;
711 if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=2999;
712 if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-1000;
715 case 200: // all light flavoured mesons (as for the 100 case)
717 result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=299;
719 if(requiredPDG>0)result = particlePDG>=100 && particlePDG<=299;
720 if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-100;
723 case 2000: // all light flavoured baryons (as for the 1000 case)
725 result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=2999;
727 if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=2999;
728 if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-1000;
731 case 300: // all strange mesons
733 result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
735 if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
736 if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
739 case 3000: // all strange baryons
741 result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
743 if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
744 if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
747 case 400: // all charmed mesons
749 result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
751 if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
752 if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
755 case 4000: // all charmed baryons
757 result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
759 if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
760 if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
763 case 500: // all beauty mesons
765 result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
767 if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
768 if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
771 case 5000: // all beauty baryons
773 result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
775 if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
776 if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
779 default: // all specific cases
781 result = (absRequiredPDG==TMath::Abs(particlePDG));
783 result = (requiredPDG==particlePDG);
790 //________________________________________________________________________________
791 Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
793 // Check the source for the particle
797 case AliDielectronSignalMC::kDontCare :
800 case AliDielectronSignalMC::kPrimary :
801 if(label>=0 && label<GetNPrimary()) return kTRUE;
804 case AliDielectronSignalMC::kSecondary :
805 if(label>=GetNPrimary()) return kTRUE;
808 case AliDielectronSignalMC::kDirect :
809 if(label>=0 && GetMothersLabel(label)<0) return kTRUE;
812 case AliDielectronSignalMC::kDecayProduct :
813 if(label>=0 && GetMothersLabel(label)>=0) return kTRUE;
823 //________________________________________________________________________________
824 Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) {
826 // Check if the particle corresponds to the MC truth in signalMC in the branch specified
829 // NOTE: Some particles have the sign of the label flipped. It is related to the quality of matching
830 // between the ESD and the MC track. The negative labels indicate a poor matching quality
831 //if(label<0) return kFALSE;
832 if(label<0) label *= -1;
833 AliVParticle* part = GetMCTrackFromMCEvent(label);
835 AliError(Form("Could not find MC particle with label %d",label));
839 if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
840 if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
843 AliVParticle* mcMother=0x0;
845 if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
847 mLabel = GetMothersLabel(label);
848 mcMother = GetMCTrackFromMCEvent(mLabel);
850 if(!mcMother) return kFALSE;
852 if(!ComparePDG(mcMother->PdgCode(),signalMC->GetMotherPDG(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
853 if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
856 // check the grandmother
857 if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
858 AliVParticle* mcGrandMother=0x0;
861 gmLabel = GetMothersLabel(mLabel);
862 mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
864 if(!mcGrandMother) return kFALSE;
866 if(!ComparePDG(mcGrandMother->PdgCode(),signalMC->GetGrandMotherPDG(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
867 if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
874 //________________________________________________________________________________
875 Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
877 // Check if the pair corresponds to the MC truth in signalMC
881 const AliVParticle * mcD1 = pair->GetFirstDaughter();
882 const AliVParticle * mcD2 = pair->GetSecondDaughter();
883 Int_t labelD1 = (mcD1 ? mcD1->GetLabel() : -1);
884 Int_t labelD2 = (mcD2 ? mcD2->GetLabel() : -1);
885 if(labelD1<0) labelD1 *= -1;
886 if(labelD2<0) labelD2 *= -1;
889 d1Pdg=GetPdgFromLabel(labelD1);
890 d2Pdg=GetPdgFromLabel(labelD2);
893 AliVParticle* mcM1=0x0;
894 AliVParticle* mcM2=0x0;
897 AliVParticle* mcG1 = 0x0;
898 AliVParticle* mcG2 = 0x0;
900 // make direct(1-1 and 2-2) and cross(1-2 and 2-1) comparisons for the whole branch
901 Bool_t directTerm = kTRUE;
903 directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetCheckBothChargesLegs(1))
904 && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
906 directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetCheckBothChargesLegs(2))
907 && CheckParticleSource(labelD2, signalMC->GetLegSource(2));
911 if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
912 labelM1 = GetMothersLabel(labelD1);
913 if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
914 directTerm = directTerm && mcM1
915 && ComparePDG(mcM1->PdgCode(),signalMC->GetMotherPDG(1),signalMC->GetCheckBothChargesMothers(1))
916 && CheckParticleSource(labelM1, signalMC->GetMotherSource(1));
920 if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
921 labelM2 = GetMothersLabel(labelD2);
922 if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
923 directTerm = directTerm && mcM2
924 && ComparePDG(mcM2->PdgCode(),signalMC->GetMotherPDG(2),signalMC->GetCheckBothChargesMothers(2))
925 && CheckParticleSource(labelM2, signalMC->GetMotherSource(2));
930 if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
931 labelG1 = GetMothersLabel(labelM1);
932 if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
933 directTerm = directTerm && mcG1
934 && ComparePDG(mcG1->PdgCode(),signalMC->GetGrandMotherPDG(1),signalMC->GetCheckBothChargesGrandMothers(1))
935 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
939 if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
940 labelG2 = GetMothersLabel(labelM2);
941 if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
942 directTerm = directTerm && mcG2
943 && ComparePDG(mcG2->PdgCode(),signalMC->GetGrandMotherPDG(2),signalMC->GetCheckBothChargesGrandMothers(2))
944 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
948 Bool_t crossTerm = kTRUE;
950 crossTerm = crossTerm && mcD2
951 && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetCheckBothChargesLegs(1))
952 && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
954 crossTerm = crossTerm && mcD1
955 && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetCheckBothChargesLegs(2))
956 && CheckParticleSource(labelD1, signalMC->GetLegSource(2));
959 if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
960 if(!mcM2 && labelD2>-1) {
961 labelM2 = GetMothersLabel(labelD2);
962 if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
964 crossTerm = crossTerm && mcM2
965 && ComparePDG(mcM2->PdgCode(),signalMC->GetMotherPDG(1),signalMC->GetCheckBothChargesMothers(1))
966 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1));
969 if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
970 if(!mcM1 && labelD1>-1) {
971 labelM1 = GetMothersLabel(labelD1);
972 if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
974 crossTerm = crossTerm && mcM1
975 && ComparePDG(mcM1->PdgCode(),signalMC->GetMotherPDG(2),signalMC->GetCheckBothChargesMothers(2))
976 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2));
980 if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
982 labelG2 = GetMothersLabel(labelM2);
983 if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
985 crossTerm = crossTerm && mcG2
986 && ComparePDG(mcG2->PdgCode(),signalMC->GetGrandMotherPDG(1),signalMC->GetCheckBothChargesGrandMothers(1))
987 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
990 if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
992 labelG1 = GetMothersLabel(labelM1);
993 if(labelG2>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
995 crossTerm = crossTerm && mcG1
996 && ComparePDG(mcG1->PdgCode(),signalMC->GetGrandMotherPDG(2),signalMC->GetCheckBothChargesGrandMothers(2))
997 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
1000 Bool_t motherRelation = kTRUE;
1001 if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
1002 motherRelation = motherRelation && HaveSameMother(pair);
1004 if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
1005 motherRelation = motherRelation && !HaveSameMother(pair);
1008 return ((directTerm || crossTerm) && motherRelation);
1013 //____________________________________________________________
1014 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
1017 // Check whether two particles have the same mother
1020 const AliVParticle * daughter1 = pair->GetFirstDaughter();
1021 const AliVParticle * daughter2 = pair->GetSecondDaughter();
1023 AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1024 AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
1025 if (!mcDaughter1 || !mcDaughter2) return 0;
1027 Int_t labelMother1=-1;
1028 Int_t labelMother2=-1;
1030 if (mcDaughter1->IsA()==AliMCParticle::Class()){
1031 labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1032 labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
1033 } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1034 labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1035 labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
1038 Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
1043 //________________________________________________________________
1044 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
1046 // return: "0" for primary jpsi
1047 // "1" for secondary jpsi (from beauty)
1048 // "2" for background
1049 if(!HaveSameMother(pair)) return 2;
1050 AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughter())->GetLabel());
1051 Int_t labelMother=-1;
1053 if (mcDaughter1->IsA()==AliMCParticle::Class()){
1054 labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1055 } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1056 labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1059 AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
1060 if(!IsMCMotherToEE(mcMother,443)) return 2;
1061 return IsJpsiPrimary(mcMother);
1064 //______________________________________________________________
1065 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
1067 // return: "0" for primary jpsi
1068 // "1" for secondary jpsi (come from B decay)
1072 if (particle->IsA()==AliMCParticle::Class()){
1073 labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1075 particle = GetMCTrackFromMCEvent(labelMoth);
1076 pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
1077 if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1078 labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1081 else if (particle->IsA()==AliAODMCParticle::Class()){
1082 labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1084 particle = GetMCTrackFromMCEvent(labelMoth);
1085 pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
1086 if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1087 labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();