AliPhysics  8d00e07 (8d00e07)
AliAnalysisTaskSEDmesonsFilterCJ.cxx
Go to the documentation of this file.
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 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A. Grelli (Utrecht University) a.grelli@uu.nl
23 // X. Zhang (LBNL) XMZhang@lbl.gov
24 // S. Aiola (Yale University) salvatore.aiola@cern.ch
25 // S. Antônio (University of São Paulo / Utrecht University) antonio.silva@cern.ch
26 // B. Trzeciak (Utrecht University) barbara.antonia.trzeciak@cern.ch
27 //-----------------------------------------------------------------------
28 
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
31 #include <TVector3.h>
32 #include <TROOT.h>
33 #include <TH3F.h>
34 #include <TRandom3.h>
35 
36 #include "AliLog.h"
38 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliAODMCHeader.h"
40 #include "AliAODHandler.h"
41 #include "AliAnalysisManager.h"
42 #include "AliLog.h"
43 #include "AliAODVertex.h"
44 #include "AliAODRecoCascadeHF.h"
46 #include "AliESDtrack.h"
47 #include "AliAODMCParticle.h"
48 #include "AliEmcalParticle.h"
49 #include "AliParticleContainer.h"
50 #include "AliAnalysisDataSlot.h"
51 #include "AliAnalysisDataContainer.h"
52 #include "AliStack.h"
53 #include "AliAnalysisVertexingHF.h"
54 #include "AliVertexingHFUtils.h"
55 
57 
59 
60 //_______________________________________________________________________________
63  fUseMCInfo(kFALSE),
64  fBuildRMEff(kFALSE),
65  fUsePythia(kFALSE),
66  fUseRejTracks(kFALSE),
67  fTrackIneff(0),
68  fUseReco(kTRUE),
69  fCandidateType(0),
70  fCandidateName(""),
71  fPDGmother(0),
72  fNProngs(0),
73  fBranchName(""),
74  fCuts(0),
75  fMinMass(0.),
76  fMaxMass(0.),
77  fInhibitTask(kFALSE),
78  fCombineDmesons(kFALSE),
79  fMultCand(kFALSE),
80  fAnalyseCand(0),
81  fRejectQuarkNotFound(kTRUE),
82  fRejectDfromB(kTRUE),
83  fKeepOnlyDfromB(kFALSE),
84  fAodEvent(0),
85  fMCHeader(0),
86  fRan(0),
87  fArrayDStartoD0pi(0),
88  fMCarray(0),
89  fCandidateArray(0),
90  fSideBandArray(0),
91  fCombinedDmesons(0),
92  fCombinedDmesonsBkg(0),
93  fMCCombinedDmesons(0),
94  fNCand(0),
95  fNSBCand(0),
96  fHistStat(0),
97  fHistNSBCandEv(0),
98  fHistNCandEv(0),
99  fHistImpParS(0),
100  fHistImpParB(0),
101  fHistPtPion(0),
102  fHistInvMassPtD(0),
103  fHistInvMassS(0),
104  fHistInvMassB(0),
105  fHistAlphaDDS(0),
106  fHistAlphaDpisS(0),
107  fHistAlphaDpiS(0),
108  fHistAlphaDKS(0),
109  fHistAlphaDDB(0),
110  fHistAlphaDpisB(0),
111  fHistAlphaDpiB(0),
112  fHistAlphaDKB(0),
113  fHistDeltaRDDS(0),
114  fHistDeltaRDpisS(0),
115  fHistDeltaRDpiS(0),
116  fHistDeltaRDKS(0),
117  fHistDeltaRDDB(0),
118  fHistDeltaRDpisB(0),
119  fHistDeltaRDpiB(0),
120  fHistDeltaRDKB(0),
121  fHistAlphaDpiR(0),
122  fHistAlphaDKR(0),
123  fHistDeltaRDpiR(0),
124  fHistDeltaRDKR(0)
125 {
126  //
127  // Default constructor
128  //
129 
130  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
131  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
132 
133  fNeedEmcalGeom = kFALSE;
134 }
135 
136 //_______________________________________________________________________________
138  AliAnalysisTaskEmcal(name, kTRUE),
139  fUseMCInfo(kFALSE),
140  fBuildRMEff(kFALSE),
141  fUsePythia(kFALSE),
142  fUseRejTracks(kFALSE),
143  fTrackIneff(0),
144  fUseReco(kTRUE),
145  fCandidateType(candtype),
146  fCandidateName(""),
147  fPDGmother(0),
148  fNProngs(0),
149  fBranchName(""),
150  fCuts(cuts),
151  fMinMass(0.),
152  fMaxMass(0.),
153  fInhibitTask(kFALSE),
154  fCombineDmesons(kFALSE),
155  fMultCand(kFALSE),
156  fAnalyseCand(0),
157  fRejectQuarkNotFound(kTRUE),
158  fRejectDfromB(kTRUE),
159  fKeepOnlyDfromB(kFALSE),
160  fAodEvent(0),
161  fMCHeader(0),
162  fRan(0),
163  fArrayDStartoD0pi(0),
164  fMCarray(0),
165  fCandidateArray(0),
166  fSideBandArray(0),
167  fCombinedDmesons(0),
168  fCombinedDmesonsBkg(0),
169  fMCCombinedDmesons(0),
170  fNCand(0),
171  fNSBCand(0),
172  fHistStat(0),
173  fHistNSBCandEv(0),
174  fHistNCandEv(0),
175  fHistImpParS(0),
176  fHistImpParB(0),
177  fHistPtPion(0),
178  fHistInvMassPtD(0),
179  fHistInvMassS(0),
180  fHistInvMassB(0),
181  fHistAlphaDDS(0),
182  fHistAlphaDpisS(0),
183  fHistAlphaDpiS(0),
184  fHistAlphaDKS(0),
185  fHistAlphaDDB(0),
186  fHistAlphaDpisB(0),
187  fHistAlphaDpiB(0),
188  fHistAlphaDKB(0),
189  fHistDeltaRDDS(0),
190  fHistDeltaRDpisS(0),
191  fHistDeltaRDpiS(0),
192  fHistDeltaRDKS(0),
193  fHistDeltaRDDB(0),
194  fHistDeltaRDpisB(0),
195  fHistDeltaRDpiB(0),
196  fHistDeltaRDKB(0),
197  fHistAlphaDpiR(0),
198  fHistAlphaDKR(0),
199  fHistDeltaRDpiR(0),
200  fHistDeltaRDKR(0)
201 {
202  //
203  // Constructor. Initialization of Inputs and Outputs
204  //
205 
206  Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
207 
208  fNeedEmcalGeom = kFALSE;
209 
210  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
211  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
212 
213  const Int_t nptbins = fCuts->GetNPtBins();
214  Float_t defaultSigmaD013[20]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040,0.040,0.040,0.040,0.040,0.040,0.040,0.040};
215 
216  switch (fCandidateType) {
217  case 0 :
218  fCandidateName = "D0";
219  fPDGmother = 421;
220  fNProngs = 2;
221  fPDGdaughters[0] = 211; // pi
222  fPDGdaughters[1] = 321; // K
223  fPDGdaughters[2] = 0; // empty
224  fPDGdaughters[3] = 0; // empty
225  fBranchName = "D0toKpi";
226  break;
227  case 1 :
228  fCandidateName = "DStar";
229  fPDGmother = 413;
230  fNProngs = 3;
231  fPDGdaughters[1] = 211; // pi soft
232  fPDGdaughters[0] = 421; // D0
233  fPDGdaughters[2] = 211; // pi fromD0
234  fPDGdaughters[3] = 321; // K from D0
235  fBranchName = "Dstar";
236 
237  if (nptbins<20) {
238  for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
239  }
240  else {
241  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
242  }
243  break;
244  default :
245  Warning("AliAnalysisTaskSEDmesonsFilterCJ", "%d not accepted!!", fCandidateType);
246  break;
247  }
248 
251 
252  DefineOutput(1, TList::Class()); // histos
253  DefineOutput(2, AliRDHFCuts::Class()); // my cuts
254  DefineOutput(3, TClonesArray::Class()); //array of candidates
255  DefineOutput(4, TClonesArray::Class()); //array of SB candidates
256  DefineOutput(5, TClonesArray::Class()); //array of candidates and event tracks
257  DefineOutput(6, TClonesArray::Class()); //array of SB candidates and event tracks
258  DefineOutput(7, TClonesArray::Class()); //array of MC D and event MC particles
259 }
260 
261 //_______________________________________________________________________________
263 {
264  //
265  // Destructor
266  //
267 
268  Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
269 
270  if (fCuts) {
271  delete fCuts;
272  fCuts = 0;
273  }
274 
275  if (fCandidateArray) {
276  delete fCandidateArray;
277  fCandidateArray = 0;
278  }
279 
280  if (fSideBandArray) {
281  delete fSideBandArray;
282  fSideBandArray = 0;
283  }
284 
285  if (fCombinedDmesons) {
286  delete fCombinedDmesons;
287  fCombinedDmesons = 0;
288  delete fMCCombinedDmesons;
289  fMCCombinedDmesons = 0;
290  }
291 
292  if (fCombinedDmesonsBkg) {
293  delete fCombinedDmesonsBkg;
295  }
296 }
297 
298 //_______________________________________________________________________________
300 {
301  //
302  // Initialization
303  //
304 
305  Info("AnalysisTaskSEDmesonsForJetCorrelations::Init()", "Entering method");
306 
307  switch (fCandidateType) {
308  case 0:
309  {
310  AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
311  copyfCutsDzero->SetName("AnalysisCutsDzero");
312  PostData(2, copyfCutsDzero); // Post the data
313  } break;
314  case 1:
315  {
316  AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
317  copyfCutsDstar->SetName("AnalysisCutsDStar");
318  PostData(2, copyfCutsDstar); // Post the cuts
319  } break;
320  default:
321  return;
322  }
323 
324  return;
325 }
326 
327 //_______________________________________________________________________________
329 {
330  //
331  // Create output objects
332  //
333 
334  Info("UserCreateOutputObjects","CreateOutputObjects of task %s", GetName());
335 
337 
338  DefineHistoForAnalysis(); // define histograms
339 
340  if (fUseReco) {
341  if (fCandidateType == kD0toKpi){
342  fCandidateArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
343  fSideBandArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
344  }
345  else if (fCandidateType == kDstartoKpipi) {
346  fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",10);
347  fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",10);
348  }
349  else {
350  AliWarning(Form("Candidate type %d not recognized!", fCandidateType));
351  return;
352  }
353  }
354  else {
355  fCandidateArray = new TClonesArray("AliAODMCParticle",10);
356  fSideBandArray = new TClonesArray("TObject",0); // not used
357  }
358 
359  if (fCombineDmesons) {
360  fCombinedDmesons = new TClonesArray("AliEmcalParticle",50);
361  fCombinedDmesonsBkg = new TClonesArray("AliEmcalParticle",50);
362  fMCCombinedDmesons = new TClonesArray("AliAODMCParticle",50);
363  }
364  else {
365  fCombinedDmesons = new TClonesArray("TObject",0); // not used
366  fCombinedDmesonsBkg = new TClonesArray("TObject",0); // not used
367  fMCCombinedDmesons = new TClonesArray("TObject",0); // not used
368  }
369 
370  fCandidateArray->SetOwner();
371  fCandidateArray->SetName(GetOutputSlot(3)->GetContainer()->GetName());
372 
373  //this is used for the DStar side bands and MC!
374  fSideBandArray->SetOwner();
375  fSideBandArray->SetName(GetOutputSlot(4)->GetContainer()->GetName());
376 
377  fCombinedDmesons->SetOwner();
378  fCombinedDmesons->SetName(GetOutputSlot(5)->GetContainer()->GetName());
379 
380  fCombinedDmesonsBkg->SetOwner();
381  fCombinedDmesonsBkg->SetName(GetOutputSlot(6)->GetContainer()->GetName());
382 
383  fMCCombinedDmesons->SetOwner();
384  fMCCombinedDmesons->SetName(GetOutputSlot(7)->GetContainer()->GetName());
385 
386  PostData(1, fOutput);
387  PostData(3, fCandidateArray);
388  PostData(4, fSideBandArray);
389  PostData(5, fCombinedDmesons);
390  PostData(6, fCombinedDmesonsBkg);
391  PostData(7, fMCCombinedDmesons);
392 
393  Info("UserCreateOutputObjects","Data posted for task %s", GetName());
394 }
395 
396 //_______________________________________________________________________________
398 {
399  //
400  // To be executed only once, for the first event
401  //
402 
403  AliDebug(2, "Entering ExecOnce()");
404 
405  if (fInhibitTask) return;
406 
407  // Load the event
408  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
409 
410  if (fAodEvent) {
411  fArrayDStartoD0pi = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(fBranchName.Data()));
412  }
413  else {
414  if (AODEvent() && IsStandardAOD()) {
415 
416  // In case there is an AOD handler writing a standard AOD, use the AOD
417  // event in memory rather than the input (ESD) event.
418  fAodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
419 
420  // in this case the branches in the deltaAOD (AliAOD.VertexingHF.root)
421  // have to taken from the AOD event hold by the AliAODExtension
422  AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
423  if(aodHandler->GetExtensions()) {
424  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
425  AliAODEvent *aodFromExt = ext->GetAOD();
426  fArrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
427  }
428  else {
429  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
430  fInhibitTask = kTRUE;
431  return;
432  }
433  }
434  }
435 
436  if (fArrayDStartoD0pi) {
437  TString objname(fArrayDStartoD0pi->GetClass()->GetName());
438  TClass cls(objname);
439  if (!cls.InheritsFrom("AliAODRecoDecayHF2Prong")) {
440  AliError(Form("%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
441  GetName(), cls.GetName(), fArrayDStartoD0pi->GetName()));
442  fInhibitTask = kTRUE;
443  fArrayDStartoD0pi = 0;
444  return;
445  }
446  }
447  else {
448  AliError(Form("Could not find array %s, skipping the event. Task '%s' will be disabled!", fBranchName.Data(), GetName()));
449  fInhibitTask = kTRUE;
450  return;
451  }
452 
453  if (fUseMCInfo) {
454  fMCarray = dynamic_cast<TClonesArray*>(fAodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
455  if (!fMCarray) {
456  AliError(Form("MC particles not found! Task '%s' will be disabled!", GetName()));
457  fInhibitTask = kTRUE;
458  return;
459  }
460 
461  fMCHeader = (AliAODMCHeader*)fAodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
462  if(!fMCHeader) {
463  AliError(Form("MC header not found! Task '%s' will be disabled!", GetName()));
464  return;
465  }
466  }
467 
468  if (fCombineDmesons) {
472  }
473 
474  fRan = new TRandom3(0);
475 
477 }
478 
479 //_______________________________________________________________________________
481 {
482  //
483  // Analysis execution
484  //
485 
486  if (fInhibitTask) return kFALSE;
487 
488  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
489  if (matchingAODdeltaAODlevel<=0) {
490  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
491  return kFALSE;
492  }
493 
494  AliDebug(2, "Entering Run()");
495 
496  //clear the TClonesArray from the previous event
497  fCandidateArray->Clear();
498  fSideBandArray->Clear();
499  fCombinedDmesons->Clear();
500  fCombinedDmesonsBkg->Clear();
501  AliDebug(2, "TClonesArray cleared");
502 
503  fHistStat->Fill(0);
504 
505  // fix for temporary bug in ESDfilter
506  // the AODs with null vertex pointer didn't pass the PhysSel
507  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
508 
509  //Event selection
510  Bool_t iseventselected = fCuts->IsEventSelected(fAodEvent);
511  if (!iseventselected) return kFALSE;
512  fHistStat->Fill(1);
513 
514  AliDebug(2, "Event selected");
515 
516  const Int_t nD = fArrayDStartoD0pi->GetEntriesFast();
517  AliDebug(2, Form("Found %d vertices", nD));
518  if (!fUseMCInfo) fHistStat->Fill(2, nD);
519 
520 
521  Int_t pdgMeson = 413;
522  if (fCandidateType == kD0toKpi) pdgMeson = 421;
523 
524  fNCand = 0;
525  fNSBCand = 0;
526 
528 
529 
530 // for MC response matrix of efficiency studies, fMultCand option only
531 if(fUseMCInfo && fBuildRMEff){
532 
533  const Int_t nDMC = fMCarray->GetEntriesFast();
534 
535  for (Int_t iMCcharm = 0; iMCcharm < nDMC; iMCcharm++){ //loop over MC D
536 
537  AliAODMCParticle* charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(iMCcharm));
538  if(TMath::Abs(charmPart->GetPdgCode()) != pdgMeson) continue;
539  if(!charmPart) continue;
540  fHistStat->Fill(2);
541 
542  Int_t origin = CheckOrigin(charmPart, fMCarray);
543  if (origin < 0) continue;
544  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
545  fHistStat->Fill(6);
546  continue;
547  }
548  if (fRejectDfromB && origin == kFromBottom) {
549  fHistStat->Fill(7);
550  continue;
551  }
552  if (fKeepOnlyDfromB && origin != kFromBottom) {
553  fHistStat->Fill(8);
554  continue;
555  }
556 
557  Int_t decay = CheckDecayChannel(charmPart, fMCarray);
558  if( TMath::Abs(charmPart->GetPdgCode()) == 413 && decay != kDecayDStartoKpipi) continue;
559  if( TMath::Abs(charmPart->GetPdgCode()) == 421 && decay != kDecayD0toKpi) continue;
560 
561  fHistStat->Fill(3);
562 
563  if (fNCand==fAnalyseCand){
564 
565  new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
566 
567  AliAODRecoDecayHF2Prong* charmCand = 0;
568  AliAODRecoCascadeHF* dstar = 0;
569 
570  // loop over reco D candidates to find a match to MC
571  Int_t isRecoD = kFALSE;
572  for (Int_t icharm = 0; icharm < nD; icharm++) {
573 
574  charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
575  if (!charmCand) continue;
576  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
577 
578  Int_t nprongs = charmCand->GetNProngs();
579  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
580 
581  if (fCandidateType == kDstartoKpipi) {
582  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
583  if (!dstar) {
584  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
585  continue;
586  }
587  }
588 
589  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
590  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
591 
592  Int_t mcLabel = 0;
593  if (fCandidateType == kDstartoKpipi) mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
594  else mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
595 
596  if(mcLabel == iMCcharm) { isRecoD = kTRUE; break; }
597  }
598 
599  if (!isRecoD) break;
600  if (!charmCand) break;
601  if (fCandidateType == kDstartoKpipi && !dstar) break;
602 
603  fHistStat->Fill(4);
604 
605  // region of interest + cuts
606  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) break;
607 
608  //candidate selected by cuts and PID
609  Int_t isSelected = 0;
610  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
611  if (!isSelected) break;
612 
613  fHistStat->Fill(5);
614 
615  if (fCandidateType == kDstartoKpipi) {
616  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
617  fHistInvMassS->Fill(dstar->DeltaInvMass());
618  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
619  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
620  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
621  }
622  else {
623  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
624  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
625  fHistInvMassS->Fill(charmCand->InvMassD0());
626  fHistInvMassS->Fill(charmCand->InvMassD0bar());
627  }
628 
629  if (fCandidateType == kDstartoKpipi) {
630  Int_t isDstar = charmPart == NULL ? 0 : 1;
631  //fill histograms of kinematics, using MC truth
632  FillDStarMCTruthKinHistos(dstar, isSelected, isDstar);
633 
634  new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
635  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
636 
637  }
638  else {
639 
640  Int_t isD0 = 0;
641  Int_t pdgCode = charmPart->GetPdgCode();
642  if (pdgCode == 421) { isD0 = 1; }
643  else if (pdgCode == -421) { isD0 = -1; }
644  //fill histograms of kinematics, using MC truth
645  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
646 
647  new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
648  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
649  }
650 
651  break;
652  }
653  else { fNCand++; }
654 
655 }
656 }
657 // for data, or MC without RM or efficiency studies
658 else {
659  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
660  Int_t isSelected = 0;
661 
662  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
663  if (!charmCand) continue;
664  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
665 
666  Int_t nprongs = charmCand->GetNProngs();
667  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
668 
669  AliAODRecoCascadeHF* dstar = 0;
670 
671  if (fCandidateType == kDstartoKpipi) {
672  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
673  if (!dstar) {
674  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
675  continue;
676  }
677  }
678 
679  // region of interest + cuts
680  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) continue;
681 
683  FillDstarSideBands(dstar);
684  }
685 
686 
687  //candidate selected by cuts and PID
688  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
689  if (!isSelected) continue;
690 
691  if (fCandidateType == kDstartoKpipi) {
692  ProcessDstar(dstar, isSelected);
693  }
694  else {
695  ProcessD0(charmCand, isSelected);
696  }
697  } // end of D cand loop
698 }
699 
700 
701  delete vHF;
702 
703 
704  AliDebug(2, "Loop done");
705 
706  if (fCombineDmesons) {
707  if (fCombinedDmesons->GetEntriesFast() > 0) {
709  }
710 
711  if (fMCCombinedDmesons->GetEntriesFast() > 0) {
713  }
714 
715  if (fCombinedDmesonsBkg->GetEntriesFast() > 0) {
717  }
718  }
719 
720  fHistNCandEv->Fill(fCandidateArray->GetEntriesFast());
722  Int_t nsbcand = fSideBandArray->GetEntriesFast();
723  fHistStat->Fill(4, nsbcand);
724  fHistNSBCandEv->Fill(nsbcand);
725  }
726 
727  PostData(1, fOutput);
728  PostData(3, fCandidateArray);
729  PostData(4, fSideBandArray);
730  PostData(5, fCombinedDmesons);
731  PostData(6, fCombinedDmesonsBkg);
732  PostData(7, fMCCombinedDmesons);
733 
734  AliDebug(2, "Exiting method");
735 
736  return kTRUE;
737 }
738 
739 //_______________________________________________________________________________
741 {
742  // Process the D0 candidate.
743 
744  // For MC analysis look for the AOD MC particle associated with the charm candidate
745  AliAODMCParticle* charmPart = 0;
746  if (fUseMCInfo) {
747  Int_t mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
748 
749  if (mcLabel >= 0) {
750  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
751  }
752  }
753 
754  if (charmPart) {
755 
756  Int_t origin = CheckOrigin(charmPart, fMCarray);
757  if (origin < 0) return;
758 
759  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
760  fHistStat->Fill(6);
761  return;
762  }
763  if (fRejectDfromB && origin == kFromBottom) {
764  fHistStat->Fill(7);
765  return;
766  }
767  if (fKeepOnlyDfromB && origin != kFromBottom) {
768  fHistStat->Fill(8);
769  return;
770  }
771 
772  fHistStat->Fill(2);
773  }
774  else {
775  fHistStat->Fill(5);
776  }
777 
778  // For MC background fill fSideBandArray
779  if (fUseMCInfo && !charmPart) {
780  if (fUseReco) {
781  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoDecayHF2Prong(*charmCand);
782  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
783 
784  if (fCombineDmesons) {
785  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(charmCand);
786  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(charmCand);
787  }
788 
789  fHistImpParB->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
790  fHistImpParB->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
791  fHistInvMassB->Fill(charmCand->InvMassD0());
792  fHistInvMassB->Fill(charmCand->InvMassD0bar());
793 
794  fNSBCand++;
795  }
796  }
797  // For data or MC signal fill fCandidateArray
798  else {
799  // For data or MC with the requirement fUseReco fill with candidates
800  if (fUseReco) {
801  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoDecayHF2Prong(*charmCand);
802  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
803 
804  if (fCombineDmesons) {
805  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(charmCand);
806  else if(fNCand==fAnalyseCand)
807  {
808  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
809  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
810  }
811  }
812 
813  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
814  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
815  fHistInvMassS->Fill(charmCand->InvMassD0());
816  fHistInvMassS->Fill(charmCand->InvMassD0bar());
817  }
818  // For MC with requirement particle level fill with AliAODMCParticle
819  else {
820  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
821  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
822  }
823  fHistStat->Fill(3);
824  fNCand++;
825  }
826 
827  // Now filling some more histograms
828 
829  // mass vs pt
830  if (isSelected == 1 || isSelected == 3) fHistInvMassPtD->Fill(charmCand->InvMassD0(), charmCand->Pt());
831  if (isSelected >= 2) fHistInvMassPtD->Fill(charmCand->InvMassD0bar(), charmCand->Pt());
832 
833  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
834  Int_t isD0 = 0;
835  if (charmPart) {
836  Int_t pdgCode = charmPart->GetPdgCode();
837 
838  if (pdgCode == 421) {
839  isD0 = 1;
840  }
841  else if (pdgCode == -421) {
842  isD0 = -1;
843  }
844  else {
845  AliDebug(2, "Not a D0/D0bar meson!");
846  return;
847  }
848  }
849 
850  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
851  }
852 }
853 
854 //_______________________________________________________________________________
856 {
857  // Process the D* candidate.
858 
859  AliDebug(2, "Entering method");
860 
861  // For MC analysis look for the AOD MC particle associated with the charm candidate
862  AliAODMCParticle* charmPart = 0;
863  if (fUseMCInfo) {
864  //D* and D0 prongs needed to MatchToMC method
865  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
866  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
867 
868  Int_t mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
869 
870  if (mcLabel >= 0) {
871  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
872  }
873  }
874 
875  if (charmPart) {
876 
877  Int_t origin = CheckOrigin(charmPart, fMCarray);
878  if (origin < 0) return;
879 
880  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
881  fHistStat->Fill(6);
882  return;
883  }
884  if (fRejectDfromB && origin == kFromBottom) {
885  fHistStat->Fill(7);
886  return;
887  }
888  if (fKeepOnlyDfromB && origin != kFromBottom) {
889  fHistStat->Fill(8);
890  return;
891  }
892 
893  fHistStat->Fill(2);
894  }
895  else {
896  fHistStat->Fill(5);
897  }
898 
899  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
900 
901  // For MC background fill fSideBandArray
902  if (fUseMCInfo && !charmPart) {
903  if (fUseReco) {
904  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
905  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
906 
907  if (fCombineDmesons) {
908  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(dstar);
909  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(dstar);
910  }
911 
912  fHistInvMassB->Fill(dstar->DeltaInvMass());
913  fHistImpParB->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
914  fHistImpParB->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
915  fHistImpParB->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
916 
917  fNSBCand++;
918  }
919  }
920  // For data and MC signal fill fCandidateArray
921  else {
922  // For data or MC signal with the requirement fUseReco fill with candidates
923  if (fUseReco) {
924  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoCascadeHF(*dstar);
925  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
926 
927  if (fCombineDmesons) {
928  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(dstar);
929  else if(fNCand==fAnalyseCand)
930  {
931  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
932  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
933  }
934  }
935  }
936  // For MC signal with requirement particle level fill with AliAODMCParticle
937  else {
938  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
939  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
940  }
941  fNCand++;
942 
943  fHistInvMassS->Fill(dstar->DeltaInvMass());
944  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
945  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
946  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
947  fHistStat->Fill(3);
948  }
949 
950  // Now filling some more histograms.
951 
952  // select D* in the D0 window.
953  // In the cut object window is loose to allow for side bands
954 
955  // retrieve the corresponding pt bin for the candidate
956  // and set the expected D0 width (x3)
957 
958  Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
959  Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
960  Double_t invMassPDG = mPDGDstar - mPDGD0;
961 
962  Double_t ptD = dstar->Pt();
963  Int_t ptDbin = fCuts->PtBin(ptD);
964  if (ptDbin < 0 || ptDbin >= fCuts->GetNPtBins()) {
965  AliError(Form("Pt %.3f out of bounds", ptD));
966  return;
967  }
968 
969  AliDebug(1, Form("Pt bin %d and sigma D0 %.4f", ptDbin, fSigmaD0[ptDbin]));
970  //consider the D* candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
971  if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[ptDbin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[ptDbin]))) {
972  AliDebug(2, "D0 mass within 3 sigma");
973 
974  // D* delta mass
975  AliDebug(2, "Filling invariant mass vs pt histogram");
976  fHistInvMassPtD->Fill(dstar->DeltaInvMass(), ptD); // 2 D slice for pt bins
977 
978  // Soft pion pt for good candidates
979  Double_t invmassDelta = dstar->DeltaInvMass();
980  if (TMath::Abs(invmassDelta - invMassPDG) < 0.0021) {
981  AliDebug(2, "Filling pion pt histogram");
982  //softpion from D* decay
983  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
984  if (softPionTrack) fHistPtPion->Fill(softPionTrack->Pt());
985  }
986  }
987 
988  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
989  Int_t isDstar = charmPart == NULL ? 0 : 1;
990  FillDStarMCTruthKinHistos(dstar, isDstar, isSelected);
991  }
992 
993  AliDebug(2, "Exiting method");
994 }
995 
996 //_______________________________________________________________________________
998 {
999  // Fill some histogram with kinematic information of the D0 candidate.
1000 
1001  Double_t ptD = charmCand->Pt();
1002  Double_t aD = charmCand->Phi();
1003  Double_t adaugh[2] = {charmCand->PhiProng(0), charmCand->PhiProng(1)};
1004  AliAODTrack* p0 = static_cast<AliAODTrack*>(charmCand->GetDaughter(0));
1005  AliAODTrack* p1 = static_cast<AliAODTrack*>(charmCand->GetDaughter(1));
1006  Float_t dR0 = DeltaR(charmCand, p0);
1007  Float_t dR1 = DeltaR(charmCand, p1);
1008 
1009  if (isD0 == 0) { //background
1010  if (isSelected == 1 || isSelected == 3) { // selected as D0
1011  fHistAlphaDKB->Fill(aD-adaugh[0],ptD);
1012  fHistAlphaDpiB->Fill(aD-adaugh[1],ptD);
1013 
1014  fHistDeltaRDKB->Fill(dR0,ptD);
1015  fHistDeltaRDpiB->Fill(dR1,ptD);
1016  }
1017  if (isSelected >= 2) { //selected as D0bar
1018  fHistAlphaDpiB->Fill(aD-adaugh[0],ptD);
1019  fHistAlphaDKB->Fill(aD-adaugh[1],ptD);
1020 
1021  fHistDeltaRDpiB->Fill(dR0,ptD);
1022  fHistDeltaRDKB->Fill(dR1,ptD);
1023  }
1024  }
1025  else if (isD0 == 1) { //D0
1026  fHistAlphaDKS->Fill(aD-adaugh[0],ptD);
1027  fHistAlphaDpiS->Fill(aD-adaugh[1],ptD);
1028 
1029  fHistDeltaRDKS->Fill(dR0,ptD);
1030  fHistDeltaRDpiS->Fill(dR1,ptD);
1031 
1032  if (isSelected == 3) { // selected as both D0bar/D0
1033  fHistAlphaDpiR->Fill(aD-adaugh[0],ptD);
1034  fHistAlphaDKR->Fill(aD-adaugh[1],ptD);
1035 
1036  fHistDeltaRDpiR->Fill(dR0,ptD);
1037  fHistDeltaRDKR->Fill(dR1,ptD);
1038  }
1039  }
1040  else if (isD0 == -1) { //D0bar
1041  fHistAlphaDKS->Fill(aD-adaugh[1],ptD);
1042  fHistAlphaDpiS->Fill(aD-adaugh[0],ptD);
1043 
1044  fHistDeltaRDKS->Fill(dR1,ptD);
1045  fHistDeltaRDpiS->Fill(dR0,ptD);
1046 
1047  if (isSelected == 3) { // selected as both D0bar/D0
1048  fHistAlphaDpiR->Fill(aD-adaugh[1],ptD);
1049  fHistAlphaDKR->Fill(aD-adaugh[0],ptD);
1050 
1051  fHistDeltaRDpiR->Fill(dR1, ptD);
1052  fHistDeltaRDKR->Fill(dR0, ptD);
1053  }
1054  }
1055 }
1056 
1057 //_______________________________________________________________________________
1059 {
1060  // Fill some histogram with kinematic information of the D0 candidate.
1061 
1062  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
1063  Double_t ptD = dstar->Pt();
1064  Double_t aD = dstar->Phi();
1065  Double_t apis= softPionTrack->Phi();
1066 
1067  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
1068  Double_t aD0 = D0fromDstar->Phi();
1069  Int_t isD0 = D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
1070  Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1);
1071  Double_t api = isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
1072  Double_t dRDD0 = DeltaR(dstar,D0fromDstar);
1073  Double_t dRDpis = DeltaR(dstar,softPionTrack);
1074  Double_t dRDpi = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)));
1075  Double_t dRDK = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)));
1076 
1077  if (isDstar) {
1078  fHistAlphaDDS ->Fill(aD-aD0, ptD);
1079  fHistAlphaDpisS->Fill(aD-apis, ptD);
1080  fHistAlphaDpiS ->Fill(aD-api, ptD);
1081  fHistAlphaDKS ->Fill(aD-aK, ptD);
1082 
1083  fHistDeltaRDDS ->Fill(dRDD0, ptD);
1084  fHistDeltaRDpisS->Fill(dRDpis, ptD);
1085  fHistDeltaRDpiS ->Fill(dRDpi, ptD);
1086  fHistDeltaRDKS ->Fill(dRDK, ptD);
1087  }
1088  else {
1089  fHistAlphaDDB ->Fill(aD-aD0, ptD);
1090  fHistAlphaDpisB->Fill(aD-apis, ptD);
1091  fHistAlphaDpiB ->Fill(aD-api, ptD);
1092  fHistAlphaDKB ->Fill(aD-aK, ptD);
1093 
1094  fHistDeltaRDDB ->Fill(dRDD0, ptD);
1095  fHistDeltaRDpisB->Fill(dRDpis, ptD);
1096  fHistDeltaRDpiB ->Fill(dRDpi, ptD);
1097  fHistDeltaRDKB ->Fill(dRDK, ptD);
1098  }
1099 }
1100 
1101 //_______________________________________________________________________________
1103 {
1104  // Fills the array of Dstar side band candidates.
1105 
1106  //select by track cuts the side band candidates (don't want mass cut)
1107  Int_t isSelected = fCuts->IsSelected(dstar, AliRDHFCuts::kTracks, fAodEvent);
1108  if (!isSelected) return;
1109 
1110  //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
1111  Int_t bin = fCuts->PtBin(dstar->Pt());
1112  if (bin < 0 || bin >= fCuts->GetNPtBins()) {
1113  AliError(Form("Pt %.3f out of bounds", dstar->Pt()));
1114  return;
1115  }
1116 
1117  const Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1118 
1119  //if data and Dstar from D0 side band
1120  if (((dstar->InvMassD0() < (mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0() > (mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/ ||
1121  ((dstar->InvMassD0() > (mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0() < (mPDGD0+10.*fSigmaD0[bin]))) /*right side band*/) {
1122 
1123  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
1124  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
1125  fNSBCand++;
1126  }
1127 }
1128 
1129 //_______________________________________________________________________________
1131 {
1132  // Set the mass limits using the PDG code and the given range.
1133 
1134  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
1135 
1136  // compute the Delta mass for the D*
1137  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
1138 
1139  fMinMass = mass - range;
1140  fMaxMass = mass + range;
1141 
1142  AliInfo(Form("Setting mass limits to %f, %f", fMinMass, fMaxMass));
1143  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) {
1144  AliError(Form("Wrong mass limits! Task '%s' will be disabled!", GetName()));
1145  fInhibitTask = kTRUE;
1146  }
1147 }
1148 
1149 //_______________________________________________________________________________
1151 {
1152  // Set the mass limits.
1153 
1154  if (uplimit>lowlimit) {
1155  fMinMass = lowlimit;
1156  fMaxMass = uplimit;
1157  }
1158  else {
1159  printf("Error! Lower limit larger than upper limit!\n");
1160  fMinMass = uplimit - uplimit*0.2;
1161  fMaxMass = uplimit;
1162  }
1163 }
1164 
1165 //_______________________________________________________________________________
1167 {
1168  // Set the D0 width for the D* analysis.
1169 
1170  if (nptbins > 30) {
1171  AliWarning("Maximum number of bins allowed is 30!");
1172  return kFALSE;
1173  }
1174 
1175  if (!width) return kFALSE;
1176  for (Int_t ipt=0; ipt < nptbins; ipt++) fSigmaD0[ipt] = width[ipt];
1177 
1178  return kTRUE;
1179 }
1180 
1181 //_______________________________________________________________________________
1183 {
1184  // Allocate the histograms for the analysis.
1185 
1186  // Statistics
1187  fHistStat = new TH1I("fHistStat", "Statistics", 9, -0.5, 8.5);
1188  fHistStat->GetXaxis()->SetBinLabel(1, "N ev anal");
1189  fHistStat->GetXaxis()->SetBinLabel(2, "N ev sel");
1190  if(fUseMCInfo && fBuildRMEff){
1191  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1192  fHistStat->GetXaxis()->SetBinLabel(4, "N Gen Sel D");
1193  fHistStat->GetXaxis()->SetBinLabel(5, "N D");
1194  fHistStat->GetXaxis()->SetBinLabel(6, "N cand sel cuts");
1195  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1196  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1197  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1198  }
1199  else {
1200  if (fUseMCInfo) {
1201  if (fUseReco) {
1202  fHistStat->GetXaxis()->SetBinLabel(3, "N D");
1203  }
1204  else {
1205  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1206  }
1207  }
1208  else {
1209  fHistStat->GetXaxis()->SetBinLabel(3, "N Cand");
1210  }
1211  fHistStat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
1212  if (fCandidateType == kDstartoKpipi) {
1213  fHistStat->GetXaxis()->SetBinLabel(5, "N side band cand");
1214  }
1215  if (fUseMCInfo) {
1216  fHistStat->GetXaxis()->SetBinLabel(6, "N Background");
1217  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1218  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1219  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1220  }
1221  }
1222 
1223  fHistStat->SetNdivisions(1);
1224  fOutput->Add(fHistStat);
1225 
1226  fHistNCandEv = new TH1F("fHistNCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0., 100.);
1227  fOutput->Add(fHistNCandEv);
1228 
1229  if ((fCandidateType == kDstartoKpipi) || fUseMCInfo) {
1230  fHistNSBCandEv = new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
1231  fOutput->Add(fHistNSBCandEv);
1232  }
1233 
1234  if (fCandidateType == kDstartoKpipi) {
1235  fHistPtPion = new TH1F("fHistPtPion", "Primary pions candidates pt", 500, 0., 10.);
1236  fHistPtPion->SetStats(kTRUE);
1237  fHistPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1238  fHistPtPion->GetYaxis()->SetTitle("entries");
1239  fOutput->Add(fHistPtPion);
1240  }
1241 
1242  // Invariant mass related histograms
1243  const Int_t nbinsmass = 200;
1244  const Int_t ptbinsD = 100;
1245  Float_t ptmin =0.;
1246  Float_t ptmax =50.;
1247  fHistInvMassPtD = new TH2F("fHistInvMassPtD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
1248  fHistInvMassPtD->SetStats(kTRUE);
1249  fHistInvMassPtD->GetXaxis()->SetTitle("mass (GeV/c)");
1250  fHistInvMassPtD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1251  fOutput->Add(fHistInvMassPtD);
1252 
1253  fHistImpParS = new TH2F("fHistImpParS", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1254  fOutput->Add(fHistImpParS);
1255 
1256  fHistInvMassS = new TH1F("fHistInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
1257  fHistInvMassS->SetStats(kTRUE);
1258  fHistInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
1259  fHistInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1260  fOutput->Add(fHistInvMassS);
1261 
1262  const Int_t nbinsalpha = 200;
1263  Float_t minalpha = -TMath::Pi();
1264  Float_t maxalpha = TMath::Pi();
1265  const Int_t nbinsdeltaR = 200;
1266  Float_t mindeltaR = 0.;
1267  Float_t maxdeltaR = 10.;
1268 
1269 
1270 
1271  if (fUseMCInfo) {
1272  fHistImpParB = new TH2F("fHistImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1273  fOutput->Add(fHistImpParB);
1274 
1275  fHistInvMassB = new TH1F("fHistInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
1276  fHistInvMassB->SetStats(kTRUE);
1277  fHistInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
1278  fHistInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1279  fOutput->Add(fHistInvMassB);
1280 
1281  if (fCandidateType == kDstartoKpipi) {
1282 
1283  fHistAlphaDDS = new TH2F("fHistAlphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1284  fHistAlphaDpisS = new TH2F("fHistAlphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1285  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1286  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1287 
1288  fHistAlphaDDB = new TH2F("fHistAlphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1289  fHistAlphaDpisB = new TH2F("fHistAlphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1290  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1291  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1292 
1293  fHistDeltaRDDS = new TH2F("fHistDeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1294  fHistDeltaRDpisS = new TH2F("fHistDeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1295  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1296  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1297 
1298  fHistDeltaRDDB = new TH2F("fHistDeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1299  fHistDeltaRDpisB = new TH2F("fHistDeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1300  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1301  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1302 
1303  fOutput->Add(fHistAlphaDDS);
1304  fOutput->Add(fHistAlphaDpisS);
1305  fOutput->Add(fHistAlphaDpiS);
1306  fOutput->Add(fHistAlphaDKS);
1307  fOutput->Add(fHistAlphaDDB);
1308  fOutput->Add(fHistAlphaDpisB);
1309  fOutput->Add(fHistAlphaDpiB);
1310  fOutput->Add(fHistAlphaDKB);
1311 
1312  fOutput->Add(fHistDeltaRDDS);
1313  fOutput->Add(fHistDeltaRDpisS);
1314  fOutput->Add(fHistDeltaRDpiS);
1315  fOutput->Add(fHistDeltaRDKS);
1316  fOutput->Add(fHistDeltaRDDB);
1317  fOutput->Add(fHistDeltaRDpisB);
1318  fOutput->Add(fHistDeltaRDpiB);
1319  fOutput->Add(fHistDeltaRDKB);
1320  }
1321  else if (fCandidateType == kD0toKpi) {
1322 
1323  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1324  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1325  fHistAlphaDpiR = new TH2F("fHistAlphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1326  fHistAlphaDKR = new TH2F("fHistAlphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1327 
1328  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1329  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1330 
1331  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1332  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1333  fHistDeltaRDpiR = new TH2F("fHistDeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1334  fHistDeltaRDKR = new TH2F("fHistDeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1335 
1336  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1337  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1338 
1339  fOutput->Add(fHistAlphaDpiS);
1340  fOutput->Add(fHistAlphaDKS);
1341  fOutput->Add(fHistAlphaDpiR);
1342  fOutput->Add(fHistAlphaDKR);
1343  fOutput->Add(fHistAlphaDpiB);
1344  fOutput->Add(fHistAlphaDKB);
1345 
1346  fOutput->Add(fHistDeltaRDpiS);
1347  fOutput->Add(fHistDeltaRDKS);
1348  fOutput->Add(fHistDeltaRDpiR);
1349  fOutput->Add(fHistDeltaRDKR);
1350  fOutput->Add(fHistDeltaRDpiB);
1351  fOutput->Add(fHistDeltaRDKB);
1352  }
1353  }
1354 
1355  return kTRUE;
1356 }
1357 
1358 //_______________________________________________________________________________
1359 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const
1360 {
1361  // Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1362 
1363  if (!p1 || !p2) return -1;
1364 
1365  Double_t phi1 = p1->Phi();
1366  Double_t eta1 = p1->Eta();
1367  Double_t phi2 = p2->Phi();
1368  Double_t eta2 = p2->Eta();
1369 
1370  Double_t dPhi = phi1 - phi2;
1371  if(dPhi < -(TMath::Pi())/2) dPhi = dPhi + TMath::TwoPi();
1372  if(dPhi > (3*(TMath::Pi()))/2) dPhi = dPhi - TMath::TwoPi();
1373 
1374  Double_t dEta = eta1 - eta2;
1375 
1376  Double_t deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1377 
1378  return deltaR;
1379 }
1380 
1381 //_______________________________________________________________________________
1383 {
1384  //
1385  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1386  //
1387 
1388  if (!tracks) return;
1389 
1390  TObjArray allDaughters(10);
1391  allDaughters.SetOwner(kFALSE);
1392 
1393  TIter next(coll);
1394  AliEmcalParticle* emcpart = 0;
1395  while ((emcpart = static_cast<AliEmcalParticle*>(next()))) {
1396  AliAODRecoDecay* reco = dynamic_cast<AliAODRecoDecay*>(emcpart->GetTrack());
1397  AliDebug(2, Form("Found a D meson candidtate with pT = %.3f, eta = %.3f, phi = %.3f", reco->Pt(), reco->Eta(), reco->Phi()));
1398  if (reco) AddDaughters(reco, allDaughters);
1399  }
1400 
1401  tracks->ResetCurrentID();
1402  AliVTrack* track = 0;
1403  Int_t n = coll->GetEntriesFast();
1404 
1405  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
1406  if(fUseMCInfo && fUsePythia){
1407  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
1408  bool isInj = IsTrackInjected(aodtrack, fMCHeader, fMCarray);
1409  if(!isInj) continue;
1410  }
1411 
1412  if (allDaughters.Remove(track) == 0) {
1413  new ((*coll)[n]) AliEmcalParticle(track);
1414  n++;
1415  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1416  }
1417  else {
1418  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1419  }
1420  }
1421 }
1422 
1423 //_______________________________________________________________________________
1425 {
1426  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1427 
1428  Int_t n = cand->GetNDaughters();
1429 
1430  //Printf("AddDaughters: the number of dauhters is %d", n);
1431 
1432  Int_t ntot = 0;
1433  Double_t pt = 0;
1434  for (Int_t i = 0; i < n; i++) {
1435  AliVTrack* track = dynamic_cast<AliVTrack*>(cand->GetDaughter(i));
1436  if (!track) continue;
1437 
1438  AliAODRecoDecay* cand2 = dynamic_cast<AliAODRecoDecay*>(track);
1439 
1440  if (cand2) {
1441  //Printf("Daughter pT = %.3f --> ", track->Pt());
1442  pt += AddDaughters(cand2, daughters);
1443  }
1444  else {
1445  if (!track->InheritsFrom("AliAODTrack")) {
1446  Printf("Warning: One of the daughters is not of type 'AliAODTrack' nor 'AliAODRecoDecay'.");
1447  continue;
1448  }
1449  //Printf("Daughter pT = %.3f", track->Pt());
1450  daughters.AddLast(track);
1451  pt += track->Pt();
1452  ntot++;
1453  }
1454  }
1455 
1456  //Printf("Total pt of the daughters = %.3f", pt);
1457 
1458  return pt;
1459 }
1460 //_______________________________________________________________________________
1462 {
1463  //
1464  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1465  //
1466 
1467  if (!mctracks) return;
1468 
1469  TObjArray allMCDaughters(10);
1470  allMCDaughters.SetOwner(kFALSE);
1471 
1472  AliAODMCParticle* mcD = (AliAODMCParticle*)coll->At(0);
1473 
1474  AddMCDaughters(mcD,allMCDaughters,fMCarray);
1475 
1476  mctracks->ResetCurrentID();
1477  AliAODMCParticle* mcpart = 0;
1478  Int_t n = coll->GetEntriesFast();
1479  while ((mcpart = static_cast<AliAODMCParticle*>(mctracks->GetNextAcceptParticle()))) {
1480  if(TMath::Abs(mcpart->Charge())==0) continue;
1481  if(fUsePythia){
1482  bool isInj = IsMCTrackInjected(mcpart, fMCHeader, fMCarray);
1483  if(!isInj) continue;
1484  }
1485 
1486 
1487  if (allMCDaughters.Remove(mcpart) == 0) {
1488  if(fUseRejTracks){
1489  if(fRan->Rndm() < fTrackIneff) continue;
1490  }
1491  new ((*coll)[n]) AliAODMCParticle(*mcpart);
1492  n++;
1493  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1494  }
1495  else {
1496  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1497  }
1498  }
1499 }
1500 
1501 //_______________________________________________________________________________
1502 Double_t AliAnalysisTaskSEDmesonsFilterCJ::AddMCDaughters(AliAODMCParticle* mcDmeson, TObjArray& mcdaughters, TClonesArray* mcArray)
1503 {
1504  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1505 
1506  Int_t n = mcDmeson->GetNDaughters();
1507  Int_t nD0 = mcDmeson->GetDaughter(0); // get label of the first daughter
1508  //Printf("AddDaughters: the number of dauhters is %d", n);
1509  Double_t pt = 0;
1510 
1511  for (Int_t i = 0; i < n; i++) {
1512  AliAODMCParticle* DDaughter = static_cast<AliAODMCParticle*>(mcArray->At(nD0+i));
1513  if (!DDaughter) continue;
1514 
1515  if (DDaughter->GetNDaughters()>0) {
1516  //Printf("Daughter pT = %.3f --> ", track->Pt());
1517  pt += AddMCDaughters(DDaughter, mcdaughters, mcArray);
1518  mcdaughters.AddLast(DDaughter);
1519  }
1520  else {
1521  //Printf("Daughter pT = %.3f", track->Pt());
1522  mcdaughters.AddLast(DDaughter);
1523  pt += DDaughter->Pt();
1524  }
1525  }
1526 
1527  //Printf("Total pt of the daughters = %.3f", pt);
1528 
1529  return pt;
1530 }
1531 //_________________________________________________________________________________________________
1532 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, TClonesArray* mcArray)
1533 {
1534  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1535 
1536  if (!mcArray) return -1;
1537 
1538  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1539  if (labDau0 < 0) return -1;
1540 
1541  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(mcArray->At(labDau0));
1542  return CheckOrigin(part, mcArray);
1543 }
1544 
1545 //_________________________________________________________________________________________________
1547 {
1548  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1549 
1550  if (!stack) return -1;
1551 
1552  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1553  if (labDau0 < 0) return -1;
1554 
1555  return CheckOrigin(labDau0, stack);
1556 }
1557 
1558 //_________________________________________________________________________________________________
1559 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODMCParticle* part, TClonesArray* mcArray)
1560 {
1561  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1562 
1563  if (!part) return -1;
1564  if (!mcArray) return -1;
1565 
1566  Int_t pdgGranma = 0;
1567  Int_t mother = part->GetMother();
1568  Int_t istep = 0;
1569  Int_t abspdgGranma = 0;
1570  Bool_t isFromB = kFALSE;
1571  Bool_t isQuarkFound = kFALSE;
1572 
1573  while (mother >= 0) {
1574  istep++;
1575  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1576  if (mcGranma != 0) {
1577  pdgGranma = mcGranma->GetPdgCode();
1578  abspdgGranma = TMath::Abs(pdgGranma);
1579  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1580  isFromB = kTRUE;
1581  }
1582 
1583  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1584  mother = mcGranma->GetMother();
1585  }
1586  else {
1587  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1588  break;
1589  }
1590  }
1591 
1592  if (isQuarkFound) {
1593  if (isFromB) {
1594  return kFromBottom;
1595  }
1596  else {
1597  return kFromCharm;
1598  }
1599  }
1600  else {
1601  return kQuarkNotFound;
1602  }
1603 }
1604 
1605 //_________________________________________________________________________________________________
1607 {
1608  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1609 
1610  if (!stack) return -1;
1611 
1612  TParticle* part = stack->Particle(ipart);
1613  if (!part) return -1;
1614 
1615  Int_t pdgGranma = 0;
1616  Int_t mother = part->GetFirstMother();
1617  Int_t istep = 0;
1618  Int_t abspdgGranma = 0;
1619  Bool_t isFromB = kFALSE;
1620  Bool_t isQuarkFound = kFALSE;
1621 
1622  while (mother >= 0) {
1623  istep++;
1624  TParticle* mcGranma = stack->Particle(mother);
1625  if (mcGranma != 0) {
1626  pdgGranma = mcGranma->GetPdgCode();
1627  abspdgGranma = TMath::Abs(pdgGranma);
1628  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1629  isFromB = kTRUE;
1630  }
1631 
1632  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1633  mother = mcGranma->GetFirstMother();
1634  }
1635  else {
1636  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1637  break;
1638  }
1639  }
1640 
1641  if (isQuarkFound) {
1642  if (isFromB) {
1643  return kFromBottom;
1644  }
1645  else {
1646  return kFromCharm;
1647  }
1648  }
1649  else {
1650  return kQuarkNotFound;
1651  }
1652 }
1653 
1654 //_________________________________________________________________________________________________
1655 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckDecayChannel(AliAODMCParticle* part, TClonesArray* mcArray)
1656 {
1657  // Determine the decay channel
1658 
1659  if (!part) return -1;
1660  if (!mcArray) return -1;
1661 
1663 
1664  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1665 
1666  if (part->GetNDaughters() == 2) {
1667 
1668  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1669  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1670 
1671  if (!d1 || !d2) {
1672  return decay;
1673  }
1674 
1675  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1676  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1677 
1678  if (absPdgPart == 421) { // D0 -> K pi
1679 
1680  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1681  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1682  decay = kDecayD0toKpi;
1683  }
1684  }
1685 
1686  if (absPdgPart == 413) { // D* -> D0 pi
1687 
1688  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1689  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1690  if (D0decay == kDecayD0toKpi) {
1691  decay = kDecayDStartoKpipi;
1692  }
1693  }
1694 
1695  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1696  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1697  if (D0decay == kDecayD0toKpi) {
1698  decay = kDecayDStartoKpipi;
1699  }
1700  }
1701  }
1702  }
1703 
1704  return decay;
1705 }
1706 
1707 //_________________________________________________________________________________________________
1709 {
1710  // Determine the decay channel
1711 
1712  if (!stack) return -1;
1713 
1714  TParticle* part = stack->Particle(ipart);
1715 
1716  if (!part) return -1;
1717 
1719 
1720  if (part->GetNDaughters() == 2) {
1721 
1722  Int_t id1 = part->GetDaughter(0);
1723  Int_t id2 = part->GetDaughter(1);
1724 
1725  TParticle* d1 = stack->Particle(id1);
1726  TParticle* d2 = stack->Particle(id2);
1727 
1728  if (!d1 || !d2) {
1729  return decay;
1730  }
1731 
1732  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1733  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1734 
1735 
1736  if (part->GetPdgCode() == 421) { // D0 -> K pi
1737 
1738  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1739  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1740  decay = kDecayD0toKpi;
1741  }
1742  }
1743 
1744  if (part->GetPdgCode() == 413) { // D* -> D0 pi
1745 
1746  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1747  Int_t D0decay = CheckDecayChannel(id1, stack);
1748  if (D0decay == kDecayD0toKpi) {
1749  decay = kDecayDStartoKpipi;
1750  }
1751  }
1752 
1753  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1754  Int_t D0decay = CheckDecayChannel(id2, stack);
1755  if (D0decay == kDecayD0toKpi) {
1756  decay = kDecayDStartoKpipi;
1757  }
1758  }
1759  }
1760  }
1761 
1762  return decay;
1763 
1764 }
1765 
1766 
1767 //_____________________________________________________________________
1768 void AliAnalysisTaskSEDmesonsFilterCJ::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1769 
1771 
1772  Int_t lab=TMath::Abs(track->GetLabel());
1773  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1774 
1775  // Int_t countControl=0;
1776 
1777  while(nameGen.IsWhitespace()){
1778  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1779  if(!mcpart){
1780  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1781  break;
1782  }
1783  Int_t mother = mcpart->GetMother();
1784  if(mother<0){
1785  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1786  break;
1787  }
1788  lab=mother;
1789  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1790  // countControl++;
1791  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1792  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1793  // break;
1794  // }
1795  }
1796 
1797  return;
1798 }
1799 //----------------------------------------------------------------------
1800 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1802  TString nameGen;
1803  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1804 
1805  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1806 
1807  return kTRUE;
1808 }
1809 
1810 
1811 //_____________________________________________________________________
1812 void AliAnalysisTaskSEDmesonsFilterCJ::GetMCTrackPrimaryGenerator(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1813 
1815 
1816  Int_t lab=TMath::Abs(track->GetLabel());
1817  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1818 
1819  // Int_t countControl=0;
1820 
1821  while(nameGen.IsWhitespace()){
1822  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1823  if(!mcpart){
1824  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1825  break;
1826  }
1827  Int_t mother = mcpart->GetMother();
1828  if(mother<0){
1829  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1830  break;
1831  }
1832  lab=mother;
1833  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1834  // countControl++;
1835  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1836  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1837  // break;
1838  // }
1839  }
1840 
1841  return;
1842 }
1843 //----------------------------------------------------------------------
1844 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsMCTrackInjected(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1846  TString nameGen;
1847  GetMCTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1848 
1849  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1850 
1851  return kTRUE;
1852 }
void ProcessDstar(AliAODRecoCascadeHF *dstar, Int_t isSelected)
Int_t pdg
static Int_t CheckOrigin(AliAODRecoDecay *cand, TClonesArray *mcArray)
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
Definition: External.C:236
TClonesArray * fCombinedDmesonsBkg
contains candidates selected by AliRDHFCuts and the rest of the event tracks
Double_t DeltaInvMass() const
void ProcessD0(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Base task in the EMCAL framework.
TClonesArray * fArrayDStartoD0pi
Random number generator.
Double_t mass
static Int_t CheckMatchingAODdeltaAODevents()
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
void FillD0MCTruthKinHistos(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected, Int_t isD0)
static Double_t AddDaughters(AliAODRecoDecay *cand, TObjArray &daughters)
void AddEventTracks(TClonesArray *coll, AliParticleContainer *tracks)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
AliVTrack * GetTrack() const
AliStack * stack
void GetMCTrackPrimaryGenerator(AliAODMCParticle *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
void FillDStarMCTruthKinHistos(AliAODRecoCascadeHF *dstar, Int_t isSelected, Int_t isDstar)
static Double_t AddMCDaughters(AliAODMCParticle *mcDmeson, TObjArray &mcdaughters, TClonesArray *mcArray)
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
Container for particles within the EMCAL framework.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Int_t fNCand
contains MC D0 and MC event particles
Bool_t IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC)
void FillDstarSideBands(AliAODRecoCascadeHF *dstar)
TClonesArray * fSideBandArray
contains candidates selected by AliRDHFCuts
int Int_t
Definition: External.C:63
Definition: External.C:204
static Int_t CheckDecayChannel(AliAODMCParticle *part, TClonesArray *mcArray)
float Float_t
Definition: External.C:68
const Double_t ptmax
AliAODTrack * GetBachelor() const
TClonesArray * fMCCombinedDmesons
contains bkg candidates selected by AliRDHFCuts and the rest of the event tracks
const Double_t ptmin
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t IsEventSelected(AliVEvent *event)
decay
Definition: HFPtSpectrum.C:41
AliEmcalList * fOutput
!output list
Int_t fNSBCand
number of selected D candidates already added to fCandidateArray
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t IsMCTrackInjected(AliAODMCParticle *track, AliAODMCHeader *header, TClonesArray *arrayMC)
Double_t InvMassD0() const
TH1 * fHistStat
number of selected side-band D candidates already added to fSideBandArray
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:248
Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
Add object to event.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:310
Int_t PtBin(Double_t pt) const
void AddMCEventTracks(TClonesArray *coll, AliParticleContainer *mctracks)
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t nptbins
void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
TClonesArray * fCombinedDmesons
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...