major dielectron update (included also the data and plotting macros for paper)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronMC.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 //#####################################################
17 //#                                                   # 
18 //#              Class AliDielectronMC                #
19 //#       Cut Class for Jpsi->e+e- analysis           #
20 //#                                                   #
21 //#   by WooJin J. Park, GSI / W.J.Park@gsi.de        #
22 //#                                                   #
23 //#####################################################
24
25 #include <AliAnalysisManager.h>
26 #include <AliAODHandler.h>
27 #include <AliESDInputHandler.h>
28 #include <AliAODInputHandler.h>
29 #include <AliMCEventHandler.h>
30 #include <AliMCEvent.h>
31 #include <AliMCParticle.h>
32 #include <AliAODMCParticle.h>
33 #include <AliStack.h>
34 #include <AliESDEvent.h>
35 #include <AliAODEvent.h>
36 #include <AliESDtrack.h>
37 #include <AliAODTrack.h>
38 #include <AliLog.h>
39
40 #include <TClonesArray.h>
41
42 #include "AliDielectronSignalMC.h"
43 #include "AliDielectronMC.h"
44
45 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
46
47 //____________________________________________________________
48 AliDielectronMC* AliDielectronMC::Instance()
49 {
50   //
51   // return pointer to singleton implementation
52   //
53   if (fgInstance) return fgInstance;
54
55   AnalysisType type=kUNSET;
56   Bool_t hasMC=kFALSE;
57   if (AliAnalysisManager::GetAnalysisManager()){
58     if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
59     else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
60
61     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
62     if(type == kESD) hasMC=mcHandler!=0x0;
63     }
64  
65   fgInstance=new AliDielectronMC(type);
66  
67   fgInstance->SetHasMC(hasMC);
68    
69   return fgInstance;
70 }
71
72 //____________________________________________________________
73 AliDielectronMC::AliDielectronMC(AnalysisType type):
74   fMCEvent(0x0),
75   fStack(0x0),
76   fAnaType(type),
77   fHasMC(kTRUE),
78   fMcArray(0x0)
79 {
80   //
81   // default constructor
82   //
83 }
84
85
86 //____________________________________________________________
87 AliDielectronMC::~AliDielectronMC()
88 {
89   //
90   // default destructor
91   //
92   
93 }
94
95 //____________________________________________________________
96 void AliDielectronMC::Initialize()
97 {
98   //
99   // initialize MC class
100   //
101   if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
102 }
103
104 //____________________________________________________________
105 Int_t AliDielectronMC::GetNMCTracks()
106 {
107   //
108   //  return the number of generated tracks from MC event
109   //
110   if(fAnaType == kESD){
111     if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
112     return fMCEvent->GetNumberOfTracks();}
113   else if(fAnaType == kAOD){
114     if(!fMcArray) { AliError("No fMcArray"); return 0; }
115     return fMcArray->GetEntriesFast();
116   }
117   return 0;
118 }
119
120 //____________________________________________________________
121 Int_t AliDielectronMC::GetNMCTracksFromStack()
122 {
123   //
124   //  return the number of generated tracks from stack
125   //
126   if (!fStack){ AliError("No fStack"); return -999; }
127   return fStack->GetNtrack();
128 }
129
130 //____________________________________________________________
131 Int_t AliDielectronMC::GetNPrimary()
132 {
133   //
134   //  return the number of primary track from MC event
135   //
136   if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
137   return fMCEvent->GetNumberOfPrimaries();
138 }
139
140 //____________________________________________________________
141 Int_t AliDielectronMC::GetNPrimaryFromStack()
142 {
143   //
144   //  return the number of primary track from stack
145   //
146   if (!fStack){ AliError("No fStack"); return -999; }
147   return fStack->GetNprimary();
148 }
149
150 //____________________________________________________________
151 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk)
152 {
153   //
154   // return MC track directly from MC event
155   //
156   if (itrk<0) return NULL;
157   AliVParticle * track=0x0;
158   if(fAnaType == kESD){
159     if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
160     track = fMCEvent->GetTrack(itrk); //  tracks from MC event (ESD)
161   } else if(fAnaType == kAOD) {
162     if (!fMcArray){ AliError("No fMCEvent"); return NULL;}
163     track = (AliVParticle*)fMcArray->At(itrk); //  tracks from MC event (AOD)
164   }
165   return track;
166 }
167
168 //____________________________________________________________
169 Bool_t AliDielectronMC::ConnectMCEvent()
170 {
171   //
172   // connect stack object from the mc handler
173   //
174   if(fAnaType == kESD){
175     fMCEvent=0x0;
176     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
177     if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
178     if (!mcHandler->InitOk() ) return kFALSE;
179     if (!mcHandler->TreeK() )  return kFALSE;
180     if (!mcHandler->TreeTR() ) return kFALSE;
181     
182     AliMCEvent* mcEvent = mcHandler->MCEvent();
183     if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
184     fMCEvent = mcEvent;
185     
186     if (!UpdateStack()) return kFALSE;
187   }
188   else if(fAnaType == kAOD)
189   {
190     fMcArray = 0x0;
191     AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
192     fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
193     if(!fMcArray) return kFALSE;
194   }
195   return kTRUE;
196 }
197
198 //____________________________________________________________
199 Bool_t AliDielectronMC::UpdateStack()
200 {
201   //
202   // update stack with new event
203   //
204   if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
205   AliStack* stack = fMCEvent->Stack();
206   if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
207   fStack = stack;
208   return kTRUE;
209 }
210
211 //____________________________________________________________
212 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
213 {
214   //
215   // return MC track
216   //
217   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
218
219   Int_t nStack = fMCEvent->GetNumberOfTracks();
220   Int_t label = TMath::Abs(_track->GetLabel());
221   if(label>nStack)return NULL;
222
223   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
224   return mctrack;
225 }
226
227
228 //____________________________________________________________
229 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
230 {
231   //
232   // return MC track
233   //
234  if(!fMcArray) { AliError("No fMCArray"); return NULL;}
235  Int_t label = _track->GetLabel();
236  if(label < 0) return NULL;
237  AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
238  return mctrack; 
239 }
240
241 //____________________________________________________________
242 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
243 {
244   //
245   // return MC track from stack
246   //
247   Int_t label = TMath::Abs(_track->GetLabel());
248   if (!fStack) AliWarning("fStack is not available. Update stack first.");
249   TParticle* mcpart = fStack->Particle(label);
250   if (!mcpart) return NULL;
251   return mcpart;
252 }
253
254 //____________________________________________________________
255 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
256 {
257   //
258   // return MC track mother
259   //
260   AliMCParticle* mcpart = GetMCTrack(_track);
261   if (!mcpart) return NULL;
262   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
263   if (!mcmother) return NULL;
264   return mcmother;
265 }
266
267 //______________________________________________________________
268 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
269 {
270  //
271  // return MC track mother
272  //
273  AliAODMCParticle* mcpart = GetMCTrack(_track);
274  if (!mcpart) return NULL;
275  if(mcpart->GetMother() < 0) return NULL;
276  AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
277  if (!mcmother) return NULL;
278  return mcmother;
279 }
280 //____________________________________________________________
281 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
282   //
283   // return MC track mother
284   //
285   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
286   return mcmother;
287 }
288
289 //____________________________________________________________
290 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
291   //
292   // return MC track mother
293   //
294   if( _particle->GetMother() < 0) return NULL;
295   AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
296   return mcmother;
297 }
298
299 //____________________________________________________________
300 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
301 {
302   //
303   // return MC track mother from stack
304   //
305   TParticle* mcpart = GetMCTrackFromStack(_track);
306   if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL; 
307   TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
308   if (!mcmother) return NULL;
309   return mcmother;
310 }
311
312 //____________________________________________________________
313 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
314 {
315   //
316   // return PDG code of the track from the MC truth info
317   //
318   AliMCParticle* mcpart = GetMCTrack(_track);
319   if (!mcpart) return -999;
320   return mcpart->PdgCode();
321 }
322
323 //__________________________________________________________
324 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
325 {
326  //
327  // return PDG code of the track from the MC truth info
328  //
329  AliAODMCParticle* mcpart = GetMCTrack(_track);
330  if (!mcpart) return -999;
331  return mcpart->PdgCode();
332 }
333
334 //____________________________________________________________
335 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
336 {
337   // 
338   // return MC PDG code from stack
339   //
340   TParticle* mcpart = GetMCTrackFromStack(_track);
341   if (!mcpart) return -999;
342   return mcpart->GetPdgCode();
343 }
344
345 //____________________________________________________________
346 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
347 {
348   //
349   // return PDG code of the mother track from the MC truth info
350   //
351   AliMCParticle* mcmother = GetMCTrackMother(_track);
352   if (!mcmother) return -999;
353   return mcmother->PdgCode();
354 }
355
356 //________________________________________________________
357 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
358 {
359   //
360   // return PDG code of the mother track from the MC truth info
361   //
362   AliAODMCParticle* mcmother = GetMCTrackMother(_track);
363   if (!mcmother) return -999;
364   return mcmother->PdgCode();
365 }
366
367 //____________________________________________________________
368 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
369 {
370   //
371   // return PDG code of the mother track from stack
372   //
373   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
374   if (!mcmother) return -999;
375   return mcmother->GetPdgCode();
376 }
377
378 //____________________________________________________________
379 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
380 {
381   //
382   // return process number of the track
383   //
384   AliMCParticle* mcpart = GetMCTrack(_track);
385   if (!mcpart) return -999;
386   return 0;
387 }
388
389 //____________________________________________________________
390 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
391 {
392   //
393   // return process number of the track
394   //
395   TParticle* mcpart = GetMCTrackFromStack(_track);
396   if (!mcpart) return -999;
397   return mcpart->GetUniqueID();
398 }
399
400 //____________________________________________________________
401 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
402 {
403   //
404   // returns the number of daughters
405   //
406   AliMCParticle *mcmother=GetMCTrackMother(track);
407   if(!mcmother||!mcmother->Particle()) return -999;
408 //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
409   return mcmother->Particle()->GetNDaughters();
410 }
411
412 //_________________________________________________________
413 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
414 {
415   //
416   // returns the number of daughters
417   //
418   AliAODMCParticle *mcmother=GetMCTrackMother(track);
419   if(!mcmother) return -999;
420   return NumberOfDaughters(mcmother);  
421
422 }
423
424 //____________________________________________________________
425 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
426 {
427   //
428   // returns the number of daughters
429   //
430   AliMCParticle *mcmother=GetMCTrackMother(particle);
431   if(!mcmother||!mcmother->Particle()) return -999;
432   //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
433   return mcmother->Particle()->GetNDaughters();
434 }
435
436 //____________________________________________________________
437 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
438 {
439   //
440   // returns the number of daughters
441   //
442   AliAODMCParticle *mcmother=GetMCTrackMother(particle);
443   if(!mcmother) return -999;
444   return mcmother->GetNDaughters();
445 }
446
447 //____________________________________________________________
448 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
449 {
450   //
451   // return process number of the mother of the track
452   //
453   AliMCParticle* mcmother = GetMCTrackMother(_track);
454   if (!mcmother) return -999;
455   return 0;
456 }
457
458 //____________________________________________________________
459 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
460 {
461   //
462   // return process number of the mother of the track
463   //
464   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
465   if (!mcmother) return -999;
466   return mcmother->GetUniqueID();
467 }
468
469 //____________________________________________________________
470 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
471 {
472   //
473   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
474   //
475   if (fAnaType==kESD && !fMCEvent) return kFALSE;
476   if (fAnaType==kAOD && !fMcArray) return kFALSE;
477   if (!particle) return kFALSE;
478   
479   if (particle->IsA()==AliMCParticle::Class()){
480     return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
481   } else if (particle->IsA()==AliAODMCParticle::Class()){
482    return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
483   } else {
484     AliError("Unknown particle type");
485   }
486   return kFALSE;
487 }
488
489 //____________________________________________________________
490 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
491 {
492   //
493   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
494   // ESD case
495   //
496   
497   //check pdg code
498   if (particle->PdgCode()!=pdgMother) return kFALSE;
499   Int_t ifirst = particle->GetFirstDaughter();
500   Int_t ilast  = particle->GetLastDaughter();
501   
502   //check number of daughters
503   if ((ilast-ifirst)!=1) return kFALSE;
504   AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
505   AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
506
507   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
508   if (firstD->Charge()>0){
509     if (firstD->PdgCode()!=-11) return kFALSE;
510     if (secondD->PdgCode()!=11) return kFALSE;
511   }else{
512     if (firstD->PdgCode()!=11) return kFALSE;
513     if (secondD->PdgCode()!=-11) return kFALSE;
514   }
515   
516   return kTRUE;
517 }
518
519 //____________________________________________________________
520 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
521 {
522   //
523   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
524   // AOD case
525   //
526
527   if (particle->GetPdgCode()!=pdgMother) return kFALSE;
528   if (particle->GetNDaughters()!=2) return kFALSE;
529  
530   Int_t ifirst = particle->GetDaughter(0);
531   Int_t ilast  = particle->GetDaughter(1);
532   
533   //check number of daughters
534   if ((ilast-ifirst)!=1) return kFALSE;
535   
536   AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
537   AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
538    
539   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
540  
541   if (firstD->Charge()>0){
542     if (firstD->GetPdgCode()!=-11) return kFALSE;
543     if (secondD->GetPdgCode()!=11) return kFALSE;
544   }else{
545     if (firstD->GetPdgCode()!=11) return kFALSE;
546     if (secondD->GetPdgCode()!=-11) return kFALSE;
547   }
548   return kTRUE;
549 }
550
551 //____________________________________________________________
552 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
553 {
554   //
555   // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
556   //
557   if (fAnaType==kESD){
558   if (!fMCEvent) return -1;
559   return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
560   }
561   else if (fAnaType==kAOD)
562   {
563   if (!fMcArray) return -1;
564   return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
565   }  
566
567   return -1;
568 }
569
570 //____________________________________________________________
571 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
572 {
573   //
574   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
575   //    have the same mother and the mother had pdg code pdgMother
576   // ESD case
577   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
578   //
579   
580   AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
581   AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
582   
583   if (!mcPart1||!mcPart2) return -1;
584   
585   Int_t lblMother1=mcPart1->GetMother();
586   Int_t lblMother2=mcPart2->GetMother();
587   
588   AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
589   if (!mcMother1) return -1;
590   if (lblMother1!=lblMother2) return -1;
591   if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
592   if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
593   if (mcMother1->PdgCode()!=pdgMother) return -1;
594   
595   return lblMother1;
596 }
597
598 //____________________________________________________________
599 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
600 {
601   //
602   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
603   //    have the same mother and the mother had pdg code pdgMother
604   // AOD case
605   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
606   //
607   AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
608   AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
609   
610   if (!mcPart1||!mcPart2) return -1;
611   
612   Int_t lblMother1=mcPart1->GetMother();
613   Int_t lblMother2=mcPart2->GetMother();
614   
615   AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
616   
617   if (!mcMother1) return -1;
618   if (lblMother1!=lblMother2) return -1;
619   if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
620   if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
621   if (mcMother1->GetPdgCode()!=pdgMother) return -1;
622   
623   return lblMother1;
624 }
625
626 //____________________________________________________________
627 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
628 {
629   //
630   // Get First two daughters of the mother
631   //
632   Int_t lblD1=-1;
633   Int_t lblD2=-1;
634   d1=0;
635   d2=0;
636   if (fAnaType==kAOD){
637     if(!fMcArray) return;
638     const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
639     lblD1=aodMother->GetDaughter(0);
640     lblD2=aodMother->GetDaughter(1);
641     d1 = (AliVParticle*)fMcArray->At(lblD1);
642     d2 = (AliVParticle*)fMcArray->At(lblD2);
643    } else if (fAnaType==kESD){
644     if (!fMCEvent) return;
645     const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
646     lblD1=aodMother->GetFirstDaughter();
647     lblD2=aodMother->GetLastDaughter();
648     d1=fMCEvent->GetTrack(lblD1);
649     d2=fMCEvent->GetTrack(lblD2);
650    }
651 }
652
653
654 //________________________________________________________________________________
655 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) {
656   //
657   //  Get the label of the mother for particle with label daughterLabel
658   //
659   if(daughterLabel<0) return -1;
660   if (fAnaType==kAOD) {
661     if(!fMcArray) return -1;
662     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
663   } else if(fAnaType==kESD) {
664     if (!fMCEvent) return -1;
665     return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
666   }
667   return -1;
668 }
669
670
671 //________________________________________________________________________________
672 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) {
673   //
674   //  Get particle code using the label from stack
675   //
676   if(label<0) return 0;
677   if(fAnaType==kAOD) {
678     if(!fMcArray) return 0;
679     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
680   } else if(fAnaType==kESD) {
681     if (!fMCEvent) return 0;
682     return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
683   }
684   return 0;
685 }
686
687
688 //________________________________________________________________________________
689 Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t checkBothCharges) const {
690   //
691   //  Test the PDG codes of particles with the required ones
692   //
693   Bool_t result = kTRUE;
694   Int_t absRequiredPDG = TMath::Abs(requiredPDG);
695   switch(absRequiredPDG) {
696   case 0:
697     result = kTRUE;    // PDG not required (any code will do fine)
698     break;
699   case 100:     // all light flavoured mesons
700     if(checkBothCharges)
701       result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=299;
702     else {
703       if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=299;
704       if(requiredPDG<0) result = particlePDG>=-299 && particlePDG<=-100;
705     }
706     break;
707   case 1000:     // all light flavoured baryons
708     if(checkBothCharges)
709       result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=2999;
710     else {
711       if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=2999;
712       if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-1000;
713     }
714     break;
715   case 200:     // all light flavoured mesons  (as for the 100 case)
716     if(checkBothCharges)
717       result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=299;
718     else {
719       if(requiredPDG>0)result = particlePDG>=100 && particlePDG<=299;
720       if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-100;
721     }
722     break;
723   case 2000:     // all light flavoured baryons (as for the 1000 case)
724     if(checkBothCharges)
725       result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=2999;
726     else {
727       if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=2999;
728       if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-1000;
729     }
730     break;
731   case 300:     // all strange mesons
732     if(checkBothCharges)
733       result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
734     else {
735       if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
736       if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
737     }
738     break;
739   case 3000:     // all strange baryons
740     if(checkBothCharges)
741       result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
742     else {
743       if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
744       if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
745     }
746     break;
747   case 400:     // all charmed mesons
748     if(checkBothCharges)
749       result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
750     else {
751       if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
752       if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
753     }
754     break;
755   case 4000:     // all charmed baryons
756     if(checkBothCharges)
757       result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
758     else {
759       if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
760       if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
761     }
762     break;
763   case 500:      // all beauty mesons
764     if(checkBothCharges)
765       result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
766     else {
767       if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
768       if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
769     }
770     break;
771   case 5000:      // all beauty baryons
772     if(checkBothCharges)
773       result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
774     else {
775       if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
776       if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
777     }
778     break;
779   default:          // all specific cases
780     if(checkBothCharges)
781       result = (absRequiredPDG==TMath::Abs(particlePDG));
782     else
783       result = (requiredPDG==particlePDG);
784   }
785
786   return result;
787 }
788
789
790 //________________________________________________________________________________
791 Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) {
792   //
793   //  Check the source for the particle 
794   //
795
796   switch (source) {
797     case AliDielectronSignalMC::kDontCare :
798       return kTRUE;
799     break;
800     case AliDielectronSignalMC::kPrimary :
801       if(label>=0 && label<GetNPrimary()) return kTRUE;
802       else return kFALSE;
803     break;
804     case AliDielectronSignalMC::kSecondary :
805       if(label>=GetNPrimary()) return kTRUE;
806       else return kFALSE;
807     break;
808     case AliDielectronSignalMC::kDirect :
809       if(label>=0 && GetMothersLabel(label)<0) return kTRUE;
810       else return kFALSE;
811     break;
812     case AliDielectronSignalMC::kDecayProduct :
813       if(label>=0 && GetMothersLabel(label)>=0) return kTRUE;
814       else return kFALSE;
815     break;
816     default :
817       return kFALSE;
818   }
819   return kFALSE;
820 }
821
822
823 //________________________________________________________________________________
824 Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) {
825   //
826   // Check if the particle corresponds to the MC truth in signalMC in the branch specified
827   //
828   
829   // NOTE:  Some particles have the sign of the label flipped. It is related to the quality of matching
830   //        between the ESD and the MC track. The negative labels indicate a poor matching quality
831   //if(label<0) return kFALSE;
832   if(label<0) label *= -1; 
833   AliVParticle* part = GetMCTrackFromMCEvent(label);
834   // check the leg
835   if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
836   if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
837
838   // check the mother
839   AliVParticle* mcMother=0x0;
840   Int_t mLabel = -1;
841   if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
842     if(part) {
843       mLabel = GetMothersLabel(label);
844       mcMother = GetMCTrackFromMCEvent(mLabel);
845     }
846     if(!mcMother) return kFALSE;
847
848     if(!ComparePDG(mcMother->PdgCode(),signalMC->GetMotherPDG(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
849     if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
850   }
851
852   // check the grandmother
853   if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
854     AliVParticle* mcGrandMother=0x0;
855     Int_t gmLabel = -1;
856     if(mcMother) {
857       gmLabel = GetMothersLabel(mLabel);
858       mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
859     }
860     if(!mcGrandMother) return kFALSE;
861     
862     if(!ComparePDG(mcGrandMother->PdgCode(),signalMC->GetGrandMotherPDG(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
863     if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
864   }
865
866   return kTRUE;
867 }
868
869
870 //________________________________________________________________________________
871 Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, AliDielectronSignalMC* signalMC) {
872   //
873   // Check if the pair corresponds to the MC truth in signalMC 
874   //
875  
876   // legs (daughters)
877   const AliVParticle * mcD1 = pair->GetFirstDaughter();
878   const AliVParticle * mcD2 = pair->GetSecondDaughter();
879   Int_t labelD1 = (mcD1 ? mcD1->GetLabel() : -1);
880   Int_t labelD2 = (mcD2 ? mcD2->GetLabel() : -1);
881   if(labelD1<0) labelD1 *= -1;
882   if(labelD2<0) labelD2 *= -1;
883   Int_t d1Pdg = 0;
884   Int_t d2Pdg = 0;
885   d1Pdg=GetPdgFromLabel(labelD1);
886   d2Pdg=GetPdgFromLabel(labelD2);
887   
888   // mothers
889   AliVParticle* mcM1=0x0;
890   AliVParticle* mcM2=0x0;
891   
892   // grand-mothers
893   AliVParticle* mcG1 = 0x0;
894   AliVParticle* mcG2 = 0x0;
895   
896   // make direct(1-1 and 2-2) and cross(1-2 and 2-1) comparisons for the whole branch
897   Bool_t directTerm = kTRUE;
898   // daughters
899   directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetCheckBothChargesLegs(1)) 
900                && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
901     
902   directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetCheckBothChargesLegs(2))
903                && CheckParticleSource(labelD2, signalMC->GetLegSource(2));
904     
905   // mothers
906   Int_t labelM1 = -1;
907   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
908     labelM1 = GetMothersLabel(labelD1);
909     if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
910     directTerm = directTerm && mcM1 
911                  && ComparePDG(mcM1->PdgCode(),signalMC->GetMotherPDG(1),signalMC->GetCheckBothChargesMothers(1))
912                  && CheckParticleSource(labelM1, signalMC->GetMotherSource(1));
913   }
914   
915   Int_t labelM2 = -1;
916   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
917     labelM2 = GetMothersLabel(labelD2);
918     if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
919     directTerm = directTerm && mcM2 
920                  && ComparePDG(mcM2->PdgCode(),signalMC->GetMotherPDG(2),signalMC->GetCheckBothChargesMothers(2))
921                  && CheckParticleSource(labelM2, signalMC->GetMotherSource(2));
922   }
923  
924   // grand-mothers
925   Int_t labelG1 = -1;
926   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
927     labelG1 = GetMothersLabel(labelM1);
928     if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
929     directTerm = directTerm && mcG1 
930                  && ComparePDG(mcG1->PdgCode(),signalMC->GetGrandMotherPDG(1),signalMC->GetCheckBothChargesGrandMothers(1))
931                  && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
932   }
933
934   Int_t labelG2 = -1;
935   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
936     labelG2 = GetMothersLabel(labelM2);
937     if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
938     directTerm = directTerm && mcG2 
939                  && ComparePDG(mcG2->PdgCode(),signalMC->GetGrandMotherPDG(2),signalMC->GetCheckBothChargesGrandMothers(2))
940                  && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
941   }
942
943   // Cross term
944   Bool_t crossTerm = kTRUE;
945   // daughters
946   crossTerm = crossTerm && mcD2 
947               && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetCheckBothChargesLegs(1))
948               && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
949     
950   crossTerm = crossTerm && mcD1 
951               && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetCheckBothChargesLegs(2))
952               && CheckParticleSource(labelD1, signalMC->GetLegSource(2));
953   
954   // mothers  
955   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
956     if(!mcM2 && labelD2>-1) {
957       labelM2 = GetMothersLabel(labelD2);
958       if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
959     }
960     crossTerm = crossTerm && mcM2 
961                 && ComparePDG(mcM2->PdgCode(),signalMC->GetMotherPDG(1),signalMC->GetCheckBothChargesMothers(1))
962                 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1));
963   }
964   
965   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
966     if(!mcM1 && labelD1>-1) {
967       labelM1 = GetMothersLabel(labelD1);
968       if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);  
969     }
970     crossTerm = crossTerm && mcM1 
971                 && ComparePDG(mcM1->PdgCode(),signalMC->GetMotherPDG(2),signalMC->GetCheckBothChargesMothers(2))
972                 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2));
973   }
974
975   // grand-mothers
976   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
977     if(!mcG2 && mcM2) {
978       labelG2 = GetMothersLabel(labelM2);
979       if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
980     }
981     crossTerm = crossTerm && mcG2 
982                 && ComparePDG(mcG2->PdgCode(),signalMC->GetGrandMotherPDG(1),signalMC->GetCheckBothChargesGrandMothers(1))
983                 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
984   }
985
986   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
987     if(!mcG1 && mcM1) {
988       labelG1 = GetMothersLabel(labelM1);
989       if(labelG2>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
990     }
991     crossTerm = crossTerm && mcG1 
992                 && ComparePDG(mcG1->PdgCode(),signalMC->GetGrandMotherPDG(2),signalMC->GetCheckBothChargesGrandMothers(2))
993                 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
994   }
995
996   Bool_t motherRelation = kTRUE;
997   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
998     motherRelation = motherRelation && HaveSameMother(pair);
999   }
1000   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
1001     motherRelation = motherRelation && !HaveSameMother(pair);
1002   }
1003  
1004   return ((directTerm || crossTerm) && motherRelation);
1005 }
1006
1007
1008
1009 //____________________________________________________________
1010 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair)
1011 {
1012   //
1013   // Check whether two particles have the same mother
1014   //
1015
1016   const AliVParticle * daughter1 = pair->GetFirstDaughter();
1017   const AliVParticle * daughter2 = pair->GetSecondDaughter();
1018
1019   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1020   AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
1021   if (!mcDaughter1 || !mcDaughter2) return 0;
1022
1023   Int_t labelMother1=-1;
1024   Int_t labelMother2=-1;
1025
1026   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1027     labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1028     labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
1029   } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1030     labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1031     labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
1032   }
1033
1034   Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
1035
1036   return sameMother;
1037 }
1038
1039 //________________________________________________________________
1040 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
1041 {
1042  // return: "0" for primary jpsi 
1043  //         "1" for secondary jpsi (from beauty)
1044  //         "2" for background  
1045  if(!HaveSameMother(pair)) return 2;
1046  AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughter())->GetLabel());
1047  Int_t labelMother=-1;
1048
1049   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1050      labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1051      } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1052      labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1053      }
1054
1055  AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
1056  if(!IsMCMotherToEE(mcMother,443)) return 2;
1057  return IsJpsiPrimary(mcMother);
1058 }
1059
1060 //______________________________________________________________
1061 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
1062 {
1063   // return: "0" for primary jpsi
1064   //         "1" for secondary jpsi (come from B decay)
1065  Int_t labelMoth=-1;
1066  Int_t pdgCode;
1067
1068  if (particle->IsA()==AliMCParticle::Class()){
1069      labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1070      while(labelMoth>0){
1071        particle = GetMCTrackFromMCEvent(labelMoth);
1072        pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
1073        if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1074        labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1075        }
1076     }
1077  else if (particle->IsA()==AliAODMCParticle::Class()){
1078      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1079      while(labelMoth>0){
1080      particle = GetMCTrackFromMCEvent(labelMoth);
1081      pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
1082      if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1083      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1084      }
1085   }
1086   return 0;
1087 }