]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliDielectronMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug in poisition corrected
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronMC.cxx
CommitLineData
b2a297fa 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>
fb7d2d99 28#include <AliAODInputHandler.h>
b2a297fa 29#include <AliMCEventHandler.h>
30#include <AliMCEvent.h>
31#include <AliMCParticle.h>
32#include <AliAODMCParticle.h>
33#include <AliStack.h>
34#include <AliESDEvent.h>
fb7d2d99 35#include <AliAODEvent.h>
b2a297fa 36#include <AliESDtrack.h>
fb7d2d99 37#include <AliAODTrack.h>
b2a297fa 38#include <AliLog.h>
39
fb7d2d99 40#include <TClonesArray.h>
5720c765 41#include <TParticle.h>
ba15fdfb 42
43#include "AliDielectronSignalMC.h"
b2a297fa 44#include "AliDielectronMC.h"
45
46AliDielectronMC* AliDielectronMC::fgInstance=0x0;
47
48//____________________________________________________________
49AliDielectronMC* AliDielectronMC::Instance()
50{
51 //
52 // return pointer to singleton implementation
53 //
54 if (fgInstance) return fgInstance;
fb7d2d99 55
b2a297fa 56 AnalysisType type=kUNSET;
3505bfad 57 Bool_t hasMC=kFALSE;
58 if (AliAnalysisManager::GetAnalysisManager()){
59 if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
fb7d2d99 60 else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
8df8e382 61
3505bfad 62 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
fb7d2d99 63 if(type == kESD) hasMC=mcHandler!=0x0;
64 }
65
b2a297fa 66 fgInstance=new AliDielectronMC(type);
fb7d2d99 67
3505bfad 68 fgInstance->SetHasMC(hasMC);
fb7d2d99 69
b2a297fa 70 return fgInstance;
71}
72
73//____________________________________________________________
74AliDielectronMC::AliDielectronMC(AnalysisType type):
75 fMCEvent(0x0),
76 fStack(0x0),
572b0139 77 fAnaType(type),
fb7d2d99 78 fHasMC(kTRUE),
79 fMcArray(0x0)
b2a297fa 80{
81 //
82 // default constructor
83 //
84}
85
86
87//____________________________________________________________
88AliDielectronMC::~AliDielectronMC()
89{
90 //
91 // default destructor
92 //
93
94}
95
96//____________________________________________________________
97void AliDielectronMC::Initialize()
98{
99 //
100 // initialize MC class
101 //
102 if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
103}
104
105//____________________________________________________________
106Int_t AliDielectronMC::GetNMCTracks()
107{
108 //
109 // return the number of generated tracks from MC event
110 //
fb7d2d99 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;
b2a297fa 119}
120
121//____________________________________________________________
122Int_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
2a14a7b1 131//____________________________________________________________
c8f0f810 132Int_t AliDielectronMC::GetNPrimary() const
2a14a7b1 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//____________________________________________________________
142Int_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
b2a297fa 151//____________________________________________________________
c8f0f810 152AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk) const
b2a297fa 153{
154 //
155 // return MC track directly from MC event
156 //
83e37742 157 if (itrk<0) return NULL;
fb7d2d99 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) {
40875e45 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;}
fb7d2d99 165 track = (AliVParticle*)fMcArray->At(itrk); // tracks from MC event (AOD)
166 }
b2a297fa 167 return track;
168}
169
170//____________________________________________________________
171Bool_t AliDielectronMC::ConnectMCEvent()
172{
173 //
174 // connect stack object from the mc handler
175 //
fb7d2d99 176 if(fAnaType == kESD){
177 fMCEvent=0x0;
178 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
40875e45 179 if (!mcHandler){ /*AliError("Could not retrive MC event handler!");*/ return kFALSE; }
fb7d2d99 180 if (!mcHandler->InitOk() ) return kFALSE;
181 if (!mcHandler->TreeK() ) return kFALSE;
182 if (!mcHandler->TreeTR() ) return kFALSE;
183
184 AliMCEvent* mcEvent = mcHandler->MCEvent();
40875e45 185 if (!mcEvent){ /*AliError("Could not retrieve MC event!");*/ return kFALSE; }
fb7d2d99 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()));
40875e45 195 if (!fMcArray){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
196 else fHasMC=kTRUE;
fb7d2d99 197 }
b2a297fa 198 return kTRUE;
199}
200
201//____________________________________________________________
202Bool_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//____________________________________________________________
8df8e382 215AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
b2a297fa 216{
217 //
218 // return MC track
219 //
9143d69f 220 if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
221
222 Int_t nStack = fMCEvent->GetNumberOfTracks();
b2a297fa 223 Int_t label = TMath::Abs(_track->GetLabel());
9143d69f 224 if(label>nStack)return NULL;
7329335b 225 if(label<0) return NULL;
b2a297fa 226 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
227 return mctrack;
228}
229
fb7d2d99 230
231//____________________________________________________________
232AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
233{
234 //
235 // return MC track
236 //
40875e45 237 if(!fMcArray) { AliError("No fMcArray"); return NULL;}
fb7d2d99 238 Int_t label = _track->GetLabel();
40875e45 239 if(label < 0 || label > fMcArray->GetEntriesFast()) return NULL;
fb7d2d99 240 AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
241 return mctrack;
242}
243
b2a297fa 244//____________________________________________________________
8df8e382 245TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 258AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
b2a297fa 259{
260 //
261 // return MC track mother
262 //
263 AliMCParticle* mcpart = GetMCTrack(_track);
264 if (!mcpart) return NULL;
7329335b 265 if(mcpart->GetMother()<0) return NULL;
b2a297fa 266 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
267 if (!mcmother) return NULL;
268 return mcmother;
269}
270
fb7d2d99 271//______________________________________________________________
272AliAODMCParticle* 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}
b2a297fa 284//____________________________________________________________
8df8e382 285AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
286 //
287 // return MC track mother
288 //
7329335b 289 if(_particle->GetMother() < 0) return NULL;
8df8e382 290 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
291 return mcmother;
292}
293
294//____________________________________________________________
295AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
296 //
297 // return MC track mother
298 //
fb7d2d99 299 if( _particle->GetMother() < 0) return NULL;
300 AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
8df8e382 301 return mcmother;
302}
303
304//____________________________________________________________
305TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 318Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
b2a297fa 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
fb7d2d99 328//__________________________________________________________
329Int_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
b2a297fa 339//____________________________________________________________
8df8e382 340Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 351Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
b2a297fa 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
fb7d2d99 361//________________________________________________________
362Int_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
b2a297fa 372//____________________________________________________________
8df8e382 373Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 384Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 395Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 406Int_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
fb7d2d99 417//_________________________________________________________
418Int_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
8df8e382 429//____________________________________________________________
430Int_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;
fb7d2d99 437 //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
8df8e382 438 return mcmother->Particle()->GetNDaughters();
439}
440
441//____________________________________________________________
442Int_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//____________________________________________________________
453Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
8df8e382 464Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
b2a297fa 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//____________________________________________________________
475Bool_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 //
fb7d2d99 480 if (fAnaType==kESD && !fMCEvent) return kFALSE;
481 if (fAnaType==kAOD && !fMcArray) return kFALSE;
8df8e382 482 if (!particle) return kFALSE;
b2a297fa 483
484 if (particle->IsA()==AliMCParticle::Class()){
485 return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
486 } else if (particle->IsA()==AliAODMCParticle::Class()){
fb7d2d99 487 return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
b2a297fa 488 } else {
489 AliError("Unknown particle type");
490 }
491 return kFALSE;
b2a297fa 492}
493
494//____________________________________________________________
495Bool_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));
8df8e382 511
512 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa 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//____________________________________________________________
525Bool_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 //
fb7d2d99 531
b2a297fa 532 if (particle->GetPdgCode()!=pdgMother) return kFALSE;
533 if (particle->GetNDaughters()!=2) return kFALSE;
fb7d2d99 534
b2a297fa 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;
fb7d2d99 540
b2a297fa 541 AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
542 AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
fb7d2d99 543
8df8e382 544 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
fb7d2d99 545
b2a297fa 546 if (firstD->Charge()>0){
b2a297fa 547 if (firstD->GetPdgCode()!=-11) return kFALSE;
548 if (secondD->GetPdgCode()!=11) return kFALSE;
fb7d2d99 549 }else{
550 if (firstD->GetPdgCode()!=11) return kFALSE;
551 if (secondD->GetPdgCode()!=-11) return kFALSE;
b2a297fa 552 }
553 return kTRUE;
554}
555
556//____________________________________________________________
a655b716 557Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
b2a297fa 558{
559 //
560 // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
561 //
fb7d2d99 562 if (fAnaType==kESD){
a655b716 563 if (!fMCEvent) return -1;
fb7d2d99 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
a655b716 572 return -1;
b2a297fa 573}
574
575//____________________________________________________________
a655b716 576Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
b2a297fa 577{
578 //
a655b716 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
b2a297fa 581 // ESD case
8df8e382 582 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa 583 //
584
585 AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
586 AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
587
a655b716 588 if (!mcPart1||!mcPart2) return -1;
b2a297fa 589
590 Int_t lblMother1=mcPart1->GetMother();
591 Int_t lblMother2=mcPart2->GetMother();
592
593 AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
a655b716 594 if (!mcMother1) return -1;
595 if (lblMother1!=lblMother2) return -1;
596 if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
572b0139 597 if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
a655b716 598 if (mcMother1->PdgCode()!=pdgMother) return -1;
b2a297fa 599
a655b716 600 return lblMother1;
b2a297fa 601}
602
603//____________________________________________________________
a655b716 604Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
b2a297fa 605{
606 //
a655b716 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
b2a297fa 609 // AOD case
8df8e382 610 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa 611 //
612 AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
613 AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
614
a655b716 615 if (!mcPart1||!mcPart2) return -1;
b2a297fa 616
617 Int_t lblMother1=mcPart1->GetMother();
618 Int_t lblMother2=mcPart2->GetMother();
619
620 AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
621
a655b716 622 if (!mcMother1) return -1;
623 if (lblMother1!=lblMother2) return -1;
624 if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
572b0139 625 if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
a655b716 626 if (mcMother1->GetPdgCode()!=pdgMother) return -1;
b2a297fa 627
a655b716 628 return lblMother1;
b2a297fa 629}
630
6551594b 631//____________________________________________________________
632void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
633{
634 //
635 // Get First two daughters of the mother
636 //
6551594b 637 Int_t lblD1=-1;
638 Int_t lblD2=-1;
639 d1=0;
640 d2=0;
6551594b 641 if (fAnaType==kAOD){
fb7d2d99 642 if(!fMcArray) return;
6551594b 643 const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
644 lblD1=aodMother->GetDaughter(0);
645 lblD2=aodMother->GetDaughter(1);
fb7d2d99 646 d1 = (AliVParticle*)fMcArray->At(lblD1);
647 d2 = (AliVParticle*)fMcArray->At(lblD2);
648 } else if (fAnaType==kESD){
649 if (!fMCEvent) return;
6551594b 650 const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
651 lblD1=aodMother->GetFirstDaughter();
652 lblD2=aodMother->GetLastDaughter();
fb7d2d99 653 d1=fMCEvent->GetTrack(lblD1);
654 d2=fMCEvent->GetTrack(lblD2);
655 }
6551594b 656}
8df8e382 657
ba15fdfb 658
659//________________________________________________________________________________
c8f0f810 660Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
ba15fdfb 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//________________________________________________________________________________
c8f0f810 677Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
ba15fdfb 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//________________________________________________________________________________
5720c765 694Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t pdgExclusion, Bool_t checkBothCharges) const {
ba15fdfb 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);
5720c765 700
ba15fdfb 701 switch(absRequiredPDG) {
702 case 0:
703 result = kTRUE; // PDG not required (any code will do fine)
704 break;
5720c765 705 case 100: // light flavoured mesons
ba15fdfb 706 if(checkBothCharges)
5720c765 707 result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=199;
ba15fdfb 708 else {
5720c765 709 if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=199;
710 if(requiredPDG<0) result = particlePDG>=-199 && particlePDG<=-100;
ba15fdfb 711 }
712 break;
5720c765 713 case 1000: // light flavoured baryons
ba15fdfb 714 if(checkBothCharges)
5720c765 715 result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=1999;
ba15fdfb 716 else {
5720c765 717 if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=1999;
718 if(requiredPDG<0) result = particlePDG>=-1999 && particlePDG<=-1000;
ba15fdfb 719 }
720 break;
5720c765 721 case 200: // light flavoured mesons
ba15fdfb 722 if(checkBothCharges)
5720c765 723 result = TMath::Abs(particlePDG)>=200 && TMath::Abs(particlePDG)<=299;
ba15fdfb 724 else {
5720c765 725 if(requiredPDG>0)result = particlePDG>=200 && particlePDG<=299;
726 if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-200;
ba15fdfb 727 }
728 break;
5720c765 729 case 2000: // light flavoured baryons
ba15fdfb 730 if(checkBothCharges)
5720c765 731 result = TMath::Abs(particlePDG)>=2000 && TMath::Abs(particlePDG)<=2999;
ba15fdfb 732 else {
5720c765 733 if(requiredPDG>0) result = particlePDG>=2000 && particlePDG<=2999;
734 if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-2000;
ba15fdfb 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;
5720c765 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;
ba15fdfb 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;
5720c765 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;
ba15fdfb 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;
5720c765 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;
ba15fdfb 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;
5720c765 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;
ba15fdfb 861 default: // all specific cases
862 if(checkBothCharges)
863 result = (absRequiredPDG==TMath::Abs(particlePDG));
864 else
865 result = (requiredPDG==particlePDG);
866 }
867
5720c765 868 if(absRequiredPDG!=0 && pdgExclusion) result = !result;
ba15fdfb 869 return result;
870}
871
872
5720c765 873//________________________________________________________________________________
874Bool_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
ba15fdfb 898//________________________________________________________________________________
c8f0f810 899Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
ba15fdfb 900 //
901 // Check the source for the particle
902 //
903
904 switch (source) {
905 case AliDielectronSignalMC::kDontCare :
906 return kTRUE;
907 break;
5720c765 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());
ba15fdfb 914 break;
5720c765 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);
ba15fdfb 925 break;
926 case AliDielectronSignalMC::kDirect :
5720c765 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));
ba15fdfb 939 break;
940 default :
941 return kFALSE;
942 }
943 return kFALSE;
944}
945
ef37a5a8 946//________________________________________________________________________________
ee877a49 947Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
ef37a5a8 948{
949 //
950 // Check if the particle has a three body decay, one being a photon
951 //
952 if(label<0) return kFALSE;
ee877a49 953
954
ef37a5a8 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;
ef37a5a8 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}
ba15fdfb 974
975//________________________________________________________________________________
ee877a49 976Bool_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//________________________________________________________________________________
994Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
ba15fdfb 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;
5720c765 1003
ba15fdfb 1004 AliVParticle* part = GetMCTrackFromMCEvent(label);
5cc8c825 1005 if (!part) {
1006 AliError(Form("Could not find MC particle with label %d",label));
1007 return kFALSE;
1008 }
5720c765 1009
ba15fdfb 1010 // check the leg
5720c765 1011 if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
ba15fdfb 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 }
5720c765 1022 if(!mcMother && !signalMC->GetMotherPDGexclude(branch)) return kFALSE;
ba15fdfb 1023
5720c765 1024 if(!ComparePDG((mcMother ? mcMother->PdgCode() : 0),signalMC->GetMotherPDG(branch),signalMC->GetMotherPDGexclude(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
ba15fdfb 1025 if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
ef37a5a8 1026
1027 //check for radiative deday
ee877a49 1028 if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
ba15fdfb 1029 }
5720c765 1030
ba15fdfb 1031 // check the grandmother
5720c765 1032 AliVParticle* mcGrandMother=0x0;
ba15fdfb 1033 if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
ba15fdfb 1034 Int_t gmLabel = -1;
1035 if(mcMother) {
1036 gmLabel = GetMothersLabel(mLabel);
1037 mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
1038 }
5720c765 1039 if(!mcGrandMother && !signalMC->GetGrandMotherPDGexclude(branch)) return kFALSE;
ba15fdfb 1040
5720c765 1041 if(!ComparePDG((mcGrandMother ? mcGrandMother->PdgCode() : 0),signalMC->GetGrandMotherPDG(branch),signalMC->GetGrandMotherPDGexclude(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
ba15fdfb 1042 if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
1043 }
1044
1045 return kTRUE;
1046}
1047
1048
5720c765 1049
ba15fdfb 1050//________________________________________________________________________________
c8f0f810 1051Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
ba15fdfb 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
5720c765 1079 directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
ba15fdfb 1080 && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
1081
5720c765 1082 directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
ba15fdfb 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);
ee877a49 1090 directTerm = directTerm && (mcM1 || signalMC->GetMotherPDGexclude(1))
5720c765 1091 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
ee877a49 1092 && CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
1093 && CheckRadiativeDecision(labelM1,signalMC);
ba15fdfb 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);
5720c765 1100 directTerm = directTerm && (mcM2 || signalMC->GetMotherPDGexclude(2))
1101 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
ee877a49 1102 && CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
1103 && CheckRadiativeDecision(labelM2,signalMC);
ba15fdfb 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);
5720c765 1111 directTerm = directTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(1))
1112 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
ba15fdfb 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);
5720c765 1120 directTerm = directTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(2))
1121 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
ba15fdfb 1122 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
1123 }
1124
1125 // Cross term
1126 Bool_t crossTerm = kTRUE;
1127 // daughters
1128 crossTerm = crossTerm && mcD2
5720c765 1129 && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
ba15fdfb 1130 && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
1131
1132 crossTerm = crossTerm && mcD1
5720c765 1133 && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
ba15fdfb 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 }
5720c765 1142 crossTerm = crossTerm && (mcM2 || signalMC->GetMotherPDGexclude(1))
1143 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
ee877a49 1144 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
1145 && CheckRadiativeDecision(labelM2,signalMC);
ba15fdfb 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 }
5720c765 1153 crossTerm = crossTerm && (mcM1 || signalMC->GetMotherPDGexclude(2))
1154 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
ee877a49 1155 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
1156 && CheckRadiativeDecision(labelM1,signalMC);
ba15fdfb 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 }
5720c765 1165 crossTerm = crossTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(1))
1166 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
ba15fdfb 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 }
5720c765 1175 crossTerm = crossTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(2))
1176 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
ba15fdfb 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
8df8e382 1193//____________________________________________________________
c8f0f810 1194Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
8df8e382 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();
5720c765 1202 if (!daughter1 || !daughter2) return 0;
8df8e382 1203
83e37742 1204 AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1205 AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
8df8e382 1206 if (!mcDaughter1 || !mcDaughter2) return 0;
1207
83e37742 1208 Int_t labelMother1=-1;
1209 Int_t labelMother2=-1;
8df8e382 1210
83e37742 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 }
8df8e382 1218
83e37742 1219 Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
8df8e382 1220
83e37742 1221 return sameMother;
8df8e382 1222}
1223
fb7d2d99 1224//________________________________________________________________
1225Int_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;
8df8e382 1233
fb7d2d99 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//______________________________________________________________
1246Int_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}