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