AliPhysics  7273240 (7273240)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 //-----------------------------------------------------------------------
26 
27 #include <TDatabasePDG.h>
28 #include <TParticle.h>
29 #include <TVector3.h>
30 #include <TROOT.h>
31 #include <TH3F.h>
32 
33 #include "AliLog.h"
35 #include "AliRDHFCutsD0toKpi.h"
36 #include "AliAODMCHeader.h"
37 #include "AliAODHandler.h"
38 #include "AliAnalysisManager.h"
39 #include "AliLog.h"
40 #include "AliAODVertex.h"
41 #include "AliAODRecoCascadeHF.h"
43 #include "AliESDtrack.h"
44 #include "AliAODMCParticle.h"
45 #include "AliEmcalParticle.h"
46 #include "AliParticleContainer.h"
47 #include "AliAnalysisDataSlot.h"
48 #include "AliAnalysisDataContainer.h"
49 #include "AliStack.h"
50 #include "AliAnalysisVertexingHF.h"
51 
53 
55 
56 //_______________________________________________________________________________
59  fUseMCInfo(kFALSE),
60  fUseReco(kTRUE),
61  fCandidateType(0),
62  fCandidateName(""),
63  fPDGmother(0),
64  fNProngs(0),
65  fBranchName(""),
66  fCuts(0),
67  fMinMass(0.),
68  fMaxMass(0.),
69  fInhibitTask(kFALSE),
70  fCombineDmesons(kFALSE),
71  fMultCand(kFALSE),
72  fAnalyseCand(0),
73  fRejectQuarkNotFound(kTRUE),
74  fRejectDfromB(kTRUE),
75  fKeepOnlyDfromB(kFALSE),
76  fAodEvent(0),
77  fArrayDStartoD0pi(0),
78  fMCarray(0),
79  fCandidateArray(0),
80  fSideBandArray(0),
81  fCombinedDmesons(0),
82  fCombinedDmesonsBkg(0),
83  fNCand(0),
84  fNSBCand(0),
85  fHistStat(0),
86  fHistNSBCandEv(0),
87  fHistNCandEv(0),
88  fHistImpParS(0),
89  fHistImpParB(0),
90  fHistPtPion(0),
91  fHistInvMassPtD(0),
92  fHistInvMassS(0),
93  fHistInvMassB(0),
94  fHistAlphaDDS(0),
95  fHistAlphaDpisS(0),
96  fHistAlphaDpiS(0),
97  fHistAlphaDKS(0),
98  fHistAlphaDDB(0),
99  fHistAlphaDpisB(0),
100  fHistAlphaDpiB(0),
101  fHistAlphaDKB(0),
102  fHistDeltaRDDS(0),
103  fHistDeltaRDpisS(0),
104  fHistDeltaRDpiS(0),
105  fHistDeltaRDKS(0),
106  fHistDeltaRDDB(0),
107  fHistDeltaRDpisB(0),
108  fHistDeltaRDpiB(0),
109  fHistDeltaRDKB(0),
110  fHistAlphaDpiR(0),
111  fHistAlphaDKR(0),
112  fHistDeltaRDpiR(0),
113  fHistDeltaRDKR(0)
114 {
115  //
116  // Default constructor
117  //
118 
119  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
120  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
121 
122  fNeedEmcalGeom = kFALSE;
123 }
124 
125 //_______________________________________________________________________________
127  AliAnalysisTaskEmcal(name, kTRUE),
128  fUseMCInfo(kFALSE),
129  fUseReco(kTRUE),
130  fCandidateType(candtype),
131  fCandidateName(""),
132  fPDGmother(0),
133  fNProngs(0),
134  fBranchName(""),
135  fCuts(cuts),
136  fMinMass(0.),
137  fMaxMass(0.),
138  fInhibitTask(kFALSE),
139  fCombineDmesons(kFALSE),
140  fMultCand(kFALSE),
141  fAnalyseCand(0),
142  fRejectQuarkNotFound(kTRUE),
143  fRejectDfromB(kTRUE),
144  fKeepOnlyDfromB(kFALSE),
145  fAodEvent(0),
146  fArrayDStartoD0pi(0),
147  fMCarray(0),
148  fCandidateArray(0),
149  fSideBandArray(0),
150  fCombinedDmesons(0),
151  fCombinedDmesonsBkg(0),
152  fNCand(0),
153  fNSBCand(0),
154  fHistStat(0),
155  fHistNSBCandEv(0),
156  fHistNCandEv(0),
157  fHistImpParS(0),
158  fHistImpParB(0),
159  fHistPtPion(0),
160  fHistInvMassPtD(0),
161  fHistInvMassS(0),
162  fHistInvMassB(0),
163  fHistAlphaDDS(0),
164  fHistAlphaDpisS(0),
165  fHistAlphaDpiS(0),
166  fHistAlphaDKS(0),
167  fHistAlphaDDB(0),
168  fHistAlphaDpisB(0),
169  fHistAlphaDpiB(0),
170  fHistAlphaDKB(0),
171  fHistDeltaRDDS(0),
172  fHistDeltaRDpisS(0),
173  fHistDeltaRDpiS(0),
174  fHistDeltaRDKS(0),
175  fHistDeltaRDDB(0),
176  fHistDeltaRDpisB(0),
177  fHistDeltaRDpiB(0),
178  fHistDeltaRDKB(0),
179  fHistAlphaDpiR(0),
180  fHistAlphaDKR(0),
181  fHistDeltaRDpiR(0),
182  fHistDeltaRDKR(0)
183 {
184  //
185  // Constructor. Initialization of Inputs and Outputs
186  //
187 
188  Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
189 
190  fNeedEmcalGeom = kFALSE;
191 
192  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
193  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
194 
195  const Int_t nptbins = fCuts->GetNPtBins();
196  Float_t defaultSigmaD013[13] = { 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 };
197 
198  switch (fCandidateType) {
199  case 0 :
200  fCandidateName = "D0";
201  fPDGmother = 421;
202  fNProngs = 2;
203  fPDGdaughters[0] = 211; // pi
204  fPDGdaughters[1] = 321; // K
205  fPDGdaughters[2] = 0; // empty
206  fPDGdaughters[3] = 0; // empty
207  fBranchName = "D0toKpi";
208  break;
209  case 1 :
210  fCandidateName = "DStar";
211  fPDGmother = 413;
212  fNProngs = 3;
213  fPDGdaughters[1] = 211; // pi soft
214  fPDGdaughters[0] = 421; // D0
215  fPDGdaughters[2] = 211; // pi fromD0
216  fPDGdaughters[3] = 321; // K from D0
217  fBranchName = "Dstar";
218 
219  if (nptbins<=13) {
220  for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
221  }
222  else {
223  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
224  }
225  break;
226  default :
227  Warning("AliAnalysisTaskSEDmesonsFilterCJ", "%d not accepted!!", fCandidateType);
228  break;
229  }
230 
233 
234  DefineOutput(1, TList::Class()); // histos
235  DefineOutput(2, AliRDHFCuts::Class()); // my cuts
236  DefineOutput(3, TClonesArray::Class()); //array of candidates
237  DefineOutput(4, TClonesArray::Class()); //array of SB candidates
238  DefineOutput(5, TClonesArray::Class()); //array of candidates and event tracks
239  DefineOutput(6, TClonesArray::Class()); //array of SB candidates and event tracks
240 }
241 
242 //_______________________________________________________________________________
244 {
245  //
246  // Destructor
247  //
248 
249  Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
250 
251  if (fCuts) {
252  delete fCuts;
253  fCuts = 0;
254  }
255 
256  if (fCandidateArray) {
257  delete fCandidateArray;
258  fCandidateArray = 0;
259  }
260 
261  if (fSideBandArray) {
262  delete fSideBandArray;
263  fSideBandArray = 0;
264  }
265 
266  if (fCombinedDmesons) {
267  delete fCombinedDmesons;
268  fCombinedDmesons = 0;
269  }
270 
271  if (fCombinedDmesonsBkg) {
272  delete fCombinedDmesonsBkg;
274  }
275 }
276 
277 //_______________________________________________________________________________
279 {
280  //
281  // Initialization
282  //
283 
284  Info("AnalysisTaskSEDmesonsForJetCorrelations::Init()", "Entering method");
285 
286  switch (fCandidateType) {
287  case 0:
288  {
289  AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
290  copyfCutsDzero->SetName("AnalysisCutsDzero");
291  PostData(2, copyfCutsDzero); // Post the data
292  } break;
293  case 1:
294  {
295  AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
296  copyfCutsDstar->SetName("AnalysisCutsDStar");
297  PostData(2, copyfCutsDstar); // Post the cuts
298  } break;
299  default:
300  return;
301  }
302 
303  return;
304 }
305 
306 //_______________________________________________________________________________
308 {
309  //
310  // Create output objects
311  //
312 
313  Info("UserCreateOutputObjects","CreateOutputObjects of task %s", GetName());
314 
316 
317  DefineHistoForAnalysis(); // define histograms
318 
319  if (fUseReco) {
320  if (fCandidateType == kD0toKpi){
321  fCandidateArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
322  fSideBandArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
323  }
324  else if (fCandidateType == kDstartoKpipi) {
325  fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",10);
326  fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",10);
327  }
328  else {
329  AliWarning(Form("Candidate type %d not recognized!", fCandidateType));
330  return;
331  }
332  }
333  else {
334  fCandidateArray = new TClonesArray("AliAODMCParticle",10);
335  fSideBandArray = new TClonesArray("TObject",0); // not used
336  }
337 
338  if (fCombineDmesons) {
339  fCombinedDmesons = new TClonesArray("AliEmcalParticle",50);
340  fCombinedDmesonsBkg = new TClonesArray("AliEmcalParticle",50);
341  }
342  else {
343  fCombinedDmesons = new TClonesArray("TObject",0); // not used
344  fCombinedDmesonsBkg = new TClonesArray("TObject",0); // not used
345  }
346 
347  fCandidateArray->SetOwner();
348  fCandidateArray->SetName(GetOutputSlot(3)->GetContainer()->GetName());
349 
350  //this is used for the DStar side bands and MC!
351  fSideBandArray->SetOwner();
352  fSideBandArray->SetName(GetOutputSlot(4)->GetContainer()->GetName());
353 
354  fCombinedDmesons->SetOwner();
355  fCombinedDmesons->SetName(GetOutputSlot(5)->GetContainer()->GetName());
356 
357  fCombinedDmesonsBkg->SetOwner();
358  fCombinedDmesonsBkg->SetName(GetOutputSlot(6)->GetContainer()->GetName());
359 
360  PostData(1, fOutput);
361  PostData(3, fCandidateArray);
362  PostData(4, fSideBandArray);
363  PostData(5, fCombinedDmesons);
364  PostData(6, fCombinedDmesonsBkg);
365 
366  Info("UserCreateOutputObjects","Data posted for task %s", GetName());
367 }
368 
369 //_______________________________________________________________________________
371 {
372  //
373  // To be executed only once, for the first event
374  //
375 
376  AliDebug(2, "Entering ExecOnce()");
377 
378  if (fInhibitTask) return;
379 
380  // Load the event
381  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
382 
383  if (fAodEvent) {
384  fArrayDStartoD0pi = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(fBranchName.Data()));
385  }
386  else {
387  if (AODEvent() && IsStandardAOD()) {
388 
389  // In case there is an AOD handler writing a standard AOD, use the AOD
390  // event in memory rather than the input (ESD) event.
391  fAodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
392 
393  // in this case the branches in the deltaAOD (AliAOD.VertexingHF.root)
394  // have to taken from the AOD event hold by the AliAODExtension
395  AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
396  if(aodHandler->GetExtensions()) {
397  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
398  AliAODEvent *aodFromExt = ext->GetAOD();
399  fArrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
400  }
401  else {
402  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
403  fInhibitTask = kTRUE;
404  return;
405  }
406  }
407  }
408 
409  if (fArrayDStartoD0pi) {
410  TString objname(fArrayDStartoD0pi->GetClass()->GetName());
411  TClass cls(objname);
412  if (!cls.InheritsFrom("AliAODRecoDecayHF2Prong")) {
413  AliError(Form("%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
414  GetName(), cls.GetName(), fArrayDStartoD0pi->GetName()));
415  fInhibitTask = kTRUE;
416  fArrayDStartoD0pi = 0;
417  return;
418  }
419  }
420  else {
421  AliError(Form("Could not find array %s, skipping the event. Task '%s' will be disabled!", fBranchName.Data(), GetName()));
422  fInhibitTask = kTRUE;
423  return;
424  }
425 
426  if (fUseMCInfo) {
427  fMCarray = dynamic_cast<TClonesArray*>(fAodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
428  if (!fMCarray) {
429  AliError(Form("MC particles not found! Task '%s' will be disabled!", GetName()));
430  fInhibitTask = kTRUE;
431  return;
432  }
433  }
434 
435  if (fCombineDmesons) {
438  }
439 
441 }
442 
443 //_______________________________________________________________________________
445 {
446  //
447  // Analysis execution
448  //
449 
450  if (fInhibitTask) return kFALSE;
451 
452  AliDebug(2, "Entering Run()");
453 
454  //clear the TClonesArray from the previous event
455  fCandidateArray->Clear();
456  fSideBandArray->Clear();
457  fCombinedDmesons->Clear();
458  fCombinedDmesonsBkg->Clear();
459  AliDebug(2, "TClonesArray cleared");
460 
461  fHistStat->Fill(0);
462 
463  // fix for temporary bug in ESDfilter
464  // the AODs with null vertex pointer didn't pass the PhysSel
465  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
466 
467  //Event selection
468  Bool_t iseventselected = fCuts->IsEventSelected(fAodEvent);
469  if (!iseventselected) return kFALSE;
470  fHistStat->Fill(1);
471 
472  AliDebug(2, "Event selected");
473 
474  const Int_t nD = fArrayDStartoD0pi->GetEntriesFast();
475  AliDebug(2, Form("Found %d vertices", nD));
476  if (!fUseMCInfo) fHistStat->Fill(2, nD);
477 
478  Int_t pdgMeson = 413;
479  if (fCandidateType == kD0toKpi) pdgMeson = 421;
480 
481  fNCand = 0;
482  fNSBCand = 0;
483 
484  //Fill the vertex info of the candidates
486 
487  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
488  Int_t isSelected = 0;
489 
490  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
491  if (!charmCand) continue;
492  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
493 
494  Int_t nprongs = charmCand->GetNProngs();
495  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
496 
497  AliAODRecoCascadeHF* dstar = 0;
498 
499  if (fCandidateType == kDstartoKpipi) {
500  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
501  if (!dstar) {
502  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
503  continue;
504  }
505  }
506 
507  // region of interest + cuts
508  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) continue;
509 
511  FillDstarSideBands(dstar);
512  }
513 
514  //candidate selected by cuts and PID
515  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
516 
517  if (!isSelected) continue;
518 
519  if (fCandidateType == kDstartoKpipi) {
520  ProcessDstar(dstar, isSelected);
521  }
522  else {
523  ProcessD0(charmCand, isSelected);
524  }
525  } // end of D cand loop
526 
527  delete vHF;
528 
529  AliDebug(2, "Loop done");
530 
531  if (fCombineDmesons) {
532  if (fCombinedDmesons->GetEntriesFast() > 0) {
534  }
535 
536  if (fCombinedDmesonsBkg->GetEntriesFast() > 0) {
538  }
539  }
540 
541  fHistNCandEv->Fill(fCandidateArray->GetEntriesFast());
543  Int_t nsbcand = fSideBandArray->GetEntriesFast();
544  fHistStat->Fill(4, nsbcand);
545  fHistNSBCandEv->Fill(nsbcand);
546  }
547 
548  PostData(1, fOutput);
549  PostData(3, fCandidateArray);
550  PostData(4, fSideBandArray);
551  PostData(5, fCombinedDmesons);
552  PostData(6, fCombinedDmesonsBkg);
553 
554  AliDebug(2, "Exiting method");
555 
556  return kTRUE;
557 }
558 
559 //_______________________________________________________________________________
561 {
562  // Process the D0 candidate.
563 
564  // For MC analysis look for the AOD MC particle associated with the charm candidate
565  AliAODMCParticle* charmPart = 0;
566  if (fUseMCInfo) {
567  Int_t mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
568 
569  if (mcLabel >= 0) {
570  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
571  }
572  }
573 
574  if (charmPart) {
575 
576  Int_t origin = CheckOrigin(charmPart, fMCarray);
577  if (origin < 0) return;
578 
579  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
580  fHistStat->Fill(6);
581  return;
582  }
583  if (fRejectDfromB && origin == kFromBottom) {
584  fHistStat->Fill(7);
585  return;
586  }
587  if (fKeepOnlyDfromB && origin != kFromBottom) {
588  fHistStat->Fill(8);
589  return;
590  }
591 
592  fHistStat->Fill(2);
593  }
594  else {
595  fHistStat->Fill(5);
596  }
597 
598  // For MC background fill fSideBandArray
599  if (fUseMCInfo && !charmPart) {
600  if (fUseReco) {
601  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoDecayHF2Prong(*charmCand);
602  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
603 
604  if (fCombineDmesons) {
605  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(charmCand);
606  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(charmCand);
607  }
608 
609  fHistImpParB->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
610  fHistImpParB->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
611  fHistInvMassB->Fill(charmCand->InvMassD0());
612  fHistInvMassB->Fill(charmCand->InvMassD0bar());
613 
614  fNSBCand++;
615  }
616  }
617  // For data or MC signal fill fCandidateArray
618  else {
619  // For data or MC with the requirement fUseReco fill with candidates
620  if (fUseReco) {
621  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoDecayHF2Prong(*charmCand);
622  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
623 
624  if (fCombineDmesons) {
625  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(charmCand);
626  else if(fNCand==fAnalyseCand) new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
627  }
628 
629  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
630  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
631  fHistInvMassS->Fill(charmCand->InvMassD0());
632  fHistInvMassS->Fill(charmCand->InvMassD0bar());
633  }
634  // For MC with requirement particle level fill with AliAODMCParticle
635  else {
636  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
637  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
638  }
639  fHistStat->Fill(3);
640  fNCand++;
641  }
642 
643  // Now filling some more histograms
644 
645  // mass vs pt
646  if (isSelected == 1 || isSelected == 3) fHistInvMassPtD->Fill(charmCand->InvMassD0(), charmCand->Pt());
647  if (isSelected >= 2) fHistInvMassPtD->Fill(charmCand->InvMassD0bar(), charmCand->Pt());
648 
649  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
650  Int_t isD0 = 0;
651  if (charmPart) {
652  Int_t pdgCode = charmPart->GetPdgCode();
653 
654  if (pdgCode == 421) {
655  isD0 = 1;
656  }
657  else if (pdgCode == -421) {
658  isD0 = -1;
659  }
660  else {
661  AliDebug(2, "Not a D0/D0bar meson!");
662  return;
663  }
664  }
665 
666  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
667  }
668 }
669 
670 //_______________________________________________________________________________
672 {
673  // Process the D* candidate.
674 
675  AliDebug(2, "Entering method");
676 
677  // For MC analysis look for the AOD MC particle associated with the charm candidate
678  AliAODMCParticle* charmPart = 0;
679  if (fUseMCInfo) {
680  //D* and D0 prongs needed to MatchToMC method
681  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
682  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
683 
684  Int_t mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
685 
686  if (mcLabel >= 0) {
687  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
688  }
689  }
690 
691  if (charmPart) {
692 
693  Int_t origin = CheckOrigin(charmPart, fMCarray);
694  if (origin < 0) return;
695 
696  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
697  fHistStat->Fill(6);
698  return;
699  }
700  if (fRejectDfromB && origin == kFromBottom) {
701  fHistStat->Fill(7);
702  return;
703  }
704  if (fKeepOnlyDfromB && origin != kFromBottom) {
705  fHistStat->Fill(8);
706  return;
707  }
708 
709  fHistStat->Fill(2);
710  }
711  else {
712  fHistStat->Fill(5);
713  }
714 
715  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
716 
717  // For MC background fill fSideBandArray
718  if (fUseMCInfo && !charmPart) {
719  if (fUseReco) {
720  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
721  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
722 
723  if (fCombineDmesons) {
724  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(dstar);
725  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(dstar);
726  }
727 
728  fHistInvMassB->Fill(dstar->DeltaInvMass());
729  fHistImpParB->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
730  fHistImpParB->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
731  fHistImpParB->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
732 
733  fNSBCand++;
734  }
735  }
736  // For data and MC signal fill fCandidateArray
737  else {
738  // For data or MC signal with the requirement fUseReco fill with candidates
739  if (fUseReco) {
740  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoCascadeHF(*dstar);
741  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
742 
743  if (fCombineDmesons) {
744  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(dstar);
745  else if(fNCand==fAnalyseCand) new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
746  }
747  }
748  // For MC signal with requirement particle level fill with AliAODMCParticle
749  else {
750  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
751  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
752  }
753  fNCand++;
754 
755  fHistInvMassS->Fill(dstar->DeltaInvMass());
756  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
757  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
758  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
759  fHistStat->Fill(3);
760  }
761 
762  // Now filling some more histograms.
763 
764  // select D* in the D0 window.
765  // In the cut object window is loose to allow for side bands
766 
767  // retrieve the corresponding pt bin for the candidate
768  // and set the expected D0 width (x3)
769 
770  Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
771  Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
772  Double_t invMassPDG = mPDGDstar - mPDGD0;
773 
774  Double_t ptD = dstar->Pt();
775  Int_t ptDbin = fCuts->PtBin(ptD);
776  if (ptDbin < 0 || ptDbin >= fCuts->GetNPtBins()) {
777  AliError(Form("Pt %.3f out of bounds", ptD));
778  return;
779  }
780 
781  AliDebug(1, Form("Pt bin %d and sigma D0 %.4f", ptDbin, fSigmaD0[ptDbin]));
782  //consider the D* candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
783  if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[ptDbin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[ptDbin]))) {
784  AliDebug(2, "D0 mass within 3 sigma");
785 
786  // D* delta mass
787  AliDebug(2, "Filling invariant mass vs pt histogram");
788  fHistInvMassPtD->Fill(dstar->DeltaInvMass(), ptD); // 2 D slice for pt bins
789 
790  // Soft pion pt for good candidates
791  Double_t invmassDelta = dstar->DeltaInvMass();
792  if (TMath::Abs(invmassDelta - invMassPDG) < 0.0021) {
793  AliDebug(2, "Filling pion pt histogram");
794  //softpion from D* decay
795  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
796  if (softPionTrack) fHistPtPion->Fill(softPionTrack->Pt());
797  }
798  }
799 
800  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
801  Int_t isDstar = charmPart == NULL ? 0 : 1;
802  FillDStarMCTruthKinHistos(dstar, isDstar, isSelected);
803  }
804 
805  AliDebug(2, "Exiting method");
806 }
807 
808 //_______________________________________________________________________________
810 {
811  // Fill some histogram with kinematic information of the D0 candidate.
812 
813  Double_t ptD = charmCand->Pt();
814  Double_t aD = charmCand->Phi();
815  Double_t adaugh[2] = {charmCand->PhiProng(0), charmCand->PhiProng(1)};
816  AliAODTrack* p0 = static_cast<AliAODTrack*>(charmCand->GetDaughter(0));
817  AliAODTrack* p1 = static_cast<AliAODTrack*>(charmCand->GetDaughter(1));
818  Float_t dR0 = DeltaR(charmCand, p0);
819  Float_t dR1 = DeltaR(charmCand, p1);
820 
821  if (isD0 == 0) { //background
822  if (isSelected == 1 || isSelected == 3) { // selected as D0
823  fHistAlphaDKB->Fill(aD-adaugh[0],ptD);
824  fHistAlphaDpiB->Fill(aD-adaugh[1],ptD);
825 
826  fHistDeltaRDKB->Fill(dR0,ptD);
827  fHistDeltaRDpiB->Fill(dR1,ptD);
828  }
829  if (isSelected >= 2) { //selected as D0bar
830  fHistAlphaDpiB->Fill(aD-adaugh[0],ptD);
831  fHistAlphaDKB->Fill(aD-adaugh[1],ptD);
832 
833  fHistDeltaRDpiB->Fill(dR0,ptD);
834  fHistDeltaRDKB->Fill(dR1,ptD);
835  }
836  }
837  else if (isD0 == 1) { //D0
838  fHistAlphaDKS->Fill(aD-adaugh[0],ptD);
839  fHistAlphaDpiS->Fill(aD-adaugh[1],ptD);
840 
841  fHistDeltaRDKS->Fill(dR0,ptD);
842  fHistDeltaRDpiS->Fill(dR1,ptD);
843 
844  if (isSelected == 3) { // selected as both D0bar/D0
845  fHistAlphaDpiR->Fill(aD-adaugh[0],ptD);
846  fHistAlphaDKR->Fill(aD-adaugh[1],ptD);
847 
848  fHistDeltaRDpiR->Fill(dR0,ptD);
849  fHistDeltaRDKR->Fill(dR1,ptD);
850  }
851  }
852  else if (isD0 == -1) { //D0bar
853  fHistAlphaDKS->Fill(aD-adaugh[1],ptD);
854  fHistAlphaDpiS->Fill(aD-adaugh[0],ptD);
855 
856  fHistDeltaRDKS->Fill(dR1,ptD);
857  fHistDeltaRDpiS->Fill(dR0,ptD);
858 
859  if (isSelected == 3) { // selected as both D0bar/D0
860  fHistAlphaDpiR->Fill(aD-adaugh[1],ptD);
861  fHistAlphaDKR->Fill(aD-adaugh[0],ptD);
862 
863  fHistDeltaRDpiR->Fill(dR1, ptD);
864  fHistDeltaRDKR->Fill(dR0, ptD);
865  }
866  }
867 }
868 
869 //_______________________________________________________________________________
871 {
872  // Fill some histogram with kinematic information of the D0 candidate.
873 
874  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
875  Double_t ptD = dstar->Pt();
876  Double_t aD = dstar->Phi();
877  Double_t apis= softPionTrack->Phi();
878 
879  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
880  Double_t aD0 = D0fromDstar->Phi();
881  Int_t isD0 = D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
882  Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1);
883  Double_t api = isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
884  Double_t dRDD0 = DeltaR(dstar,D0fromDstar);
885  Double_t dRDpis = DeltaR(dstar,softPionTrack);
886  Double_t dRDpi = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)));
887  Double_t dRDK = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)));
888 
889  if (isDstar) {
890  fHistAlphaDDS ->Fill(aD-aD0, ptD);
891  fHistAlphaDpisS->Fill(aD-apis, ptD);
892  fHistAlphaDpiS ->Fill(aD-api, ptD);
893  fHistAlphaDKS ->Fill(aD-aK, ptD);
894 
895  fHistDeltaRDDS ->Fill(dRDD0, ptD);
896  fHistDeltaRDpisS->Fill(dRDpis, ptD);
897  fHistDeltaRDpiS ->Fill(dRDpi, ptD);
898  fHistDeltaRDKS ->Fill(dRDK, ptD);
899  }
900  else {
901  fHistAlphaDDB ->Fill(aD-aD0, ptD);
902  fHistAlphaDpisB->Fill(aD-apis, ptD);
903  fHistAlphaDpiB ->Fill(aD-api, ptD);
904  fHistAlphaDKB ->Fill(aD-aK, ptD);
905 
906  fHistDeltaRDDB ->Fill(dRDD0, ptD);
907  fHistDeltaRDpisB->Fill(dRDpis, ptD);
908  fHistDeltaRDpiB ->Fill(dRDpi, ptD);
909  fHistDeltaRDKB ->Fill(dRDK, ptD);
910  }
911 }
912 
913 //_______________________________________________________________________________
915 {
916  // Fills the array of Dstar side band candidates.
917 
918  //select by track cuts the side band candidates (don't want mass cut)
919  Int_t isSelected = fCuts->IsSelected(dstar, AliRDHFCuts::kTracks, fAodEvent);
920  if (!isSelected) return;
921 
922  //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
923  Int_t bin = fCuts->PtBin(dstar->Pt());
924  if (bin < 0 || bin >= fCuts->GetNPtBins()) {
925  AliError(Form("Pt %.3f out of bounds", dstar->Pt()));
926  return;
927  }
928 
929  const Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
930 
931  //if data and Dstar from D0 side band
932  if (((dstar->InvMassD0() < (mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0() > (mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/ ||
933  ((dstar->InvMassD0() > (mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0() < (mPDGD0+10.*fSigmaD0[bin]))) /*right side band*/) {
934 
935  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
936  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
937  fNSBCand++;
938  }
939 }
940 
941 //_______________________________________________________________________________
943 {
944  // Set the mass limits using the PDG code and the given range.
945 
946  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
947 
948  // compute the Delta mass for the D*
949  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
950 
951  fMinMass = mass - range;
952  fMaxMass = mass + range;
953 
954  AliInfo(Form("Setting mass limits to %f, %f", fMinMass, fMaxMass));
955  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) {
956  AliError(Form("Wrong mass limits! Task '%s' will be disabled!", GetName()));
957  fInhibitTask = kTRUE;
958  }
959 }
960 
961 //_______________________________________________________________________________
963 {
964  // Set the mass limits.
965 
966  if (uplimit>lowlimit) {
967  fMinMass = lowlimit;
968  fMaxMass = uplimit;
969  }
970  else {
971  printf("Error! Lower limit larger than upper limit!\n");
972  fMinMass = uplimit - uplimit*0.2;
973  fMaxMass = uplimit;
974  }
975 }
976 
977 //_______________________________________________________________________________
979 {
980  // Set the D0 width for the D* analysis.
981 
982  if (nptbins > 30) {
983  AliWarning("Maximum number of bins allowed is 30!");
984  return kFALSE;
985  }
986 
987  if (!width) return kFALSE;
988  for (Int_t ipt=0; ipt < nptbins; ipt++) fSigmaD0[ipt] = width[ipt];
989 
990  return kTRUE;
991 }
992 
993 //_______________________________________________________________________________
995 {
996  // Allocate the histograms for the analysis.
997 
998  // Statistics
999  fHistStat = new TH1I("fHistStat", "Statistics", 9, -0.5, 8.5);
1000  fHistStat->GetXaxis()->SetBinLabel(1, "N ev anal");
1001  fHistStat->GetXaxis()->SetBinLabel(2, "N ev sel");
1002  if (fUseMCInfo) {
1003  if (fUseReco) {
1004  fHistStat->GetXaxis()->SetBinLabel(3, "N D");
1005  }
1006  else {
1007  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1008  }
1009  }
1010  else {
1011  fHistStat->GetXaxis()->SetBinLabel(3, "N Cand");
1012  }
1013  fHistStat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
1014  if (fCandidateType == kDstartoKpipi) {
1015  fHistStat->GetXaxis()->SetBinLabel(5, "N side band cand");
1016  }
1017  if (fUseMCInfo) {
1018  fHistStat->GetXaxis()->SetBinLabel(6, "N Background");
1019  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1020  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1021  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1022  }
1023 
1024  fHistStat->SetNdivisions(1);
1025  fOutput->Add(fHistStat);
1026 
1027  fHistNCandEv = new TH1F("fHistNCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0., 100.);
1028  fOutput->Add(fHistNCandEv);
1029 
1030  if ((fCandidateType == kDstartoKpipi) || fUseMCInfo) {
1031  fHistNSBCandEv = new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
1032  fOutput->Add(fHistNSBCandEv);
1033  }
1034 
1035  if (fCandidateType == kDstartoKpipi) {
1036  fHistPtPion = new TH1F("fHistPtPion", "Primary pions candidates pt", 500, 0., 10.);
1037  fHistPtPion->SetStats(kTRUE);
1038  fHistPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1039  fHistPtPion->GetYaxis()->SetTitle("entries");
1040  fOutput->Add(fHistPtPion);
1041  }
1042 
1043  // Invariant mass related histograms
1044  const Int_t nbinsmass = 200;
1045  const Int_t ptbinsD = 100;
1046  Float_t ptmin =0.;
1047  Float_t ptmax =50.;
1048  fHistInvMassPtD = new TH2F("fHistInvMassPtD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
1049  fHistInvMassPtD->SetStats(kTRUE);
1050  fHistInvMassPtD->GetXaxis()->SetTitle("mass (GeV/c)");
1051  fHistInvMassPtD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1052  fOutput->Add(fHistInvMassPtD);
1053 
1054  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
1055  fOutput->Add(fHistImpParS);
1056 
1057  fHistInvMassS = new TH1F("fHistInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
1058  fHistInvMassS->SetStats(kTRUE);
1059  fHistInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
1060  fHistInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1061  fOutput->Add(fHistInvMassS);
1062 
1063  const Int_t nbinsalpha = 200;
1064  Float_t minalpha = -TMath::Pi();
1065  Float_t maxalpha = TMath::Pi();
1066  const Int_t nbinsdeltaR = 200;
1067  Float_t mindeltaR = 0.;
1068  Float_t maxdeltaR = 10.;
1069 
1070  if (fUseMCInfo) {
1071  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
1072  fOutput->Add(fHistImpParB);
1073 
1074  fHistInvMassB = new TH1F("fHistInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
1075  fHistInvMassB->SetStats(kTRUE);
1076  fHistInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
1077  fHistInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1078  fOutput->Add(fHistInvMassB);
1079 
1080  if (fCandidateType == kDstartoKpipi) {
1081 
1082  fHistAlphaDDS = new TH2F("fHistAlphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1083  fHistAlphaDpisS = new TH2F("fHistAlphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1084  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1085  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1086 
1087  fHistAlphaDDB = new TH2F("fHistAlphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1088  fHistAlphaDpisB = new TH2F("fHistAlphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1089  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1090  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1091 
1092  fHistDeltaRDDS = new TH2F("fHistDeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1093  fHistDeltaRDpisS = new TH2F("fHistDeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1094  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1095  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1096 
1097  fHistDeltaRDDB = new TH2F("fHistDeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1098  fHistDeltaRDpisB = new TH2F("fHistDeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1099  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1100  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1101 
1102  fOutput->Add(fHistAlphaDDS);
1103  fOutput->Add(fHistAlphaDpisS);
1104  fOutput->Add(fHistAlphaDpiS);
1105  fOutput->Add(fHistAlphaDKS);
1106  fOutput->Add(fHistAlphaDDB);
1107  fOutput->Add(fHistAlphaDpisB);
1108  fOutput->Add(fHistAlphaDpiB);
1109  fOutput->Add(fHistAlphaDKB);
1110 
1111  fOutput->Add(fHistDeltaRDDS);
1112  fOutput->Add(fHistDeltaRDpisS);
1113  fOutput->Add(fHistDeltaRDpiS);
1114  fOutput->Add(fHistDeltaRDKS);
1115  fOutput->Add(fHistDeltaRDDB);
1116  fOutput->Add(fHistDeltaRDpisB);
1117  fOutput->Add(fHistDeltaRDpiB);
1118  fOutput->Add(fHistDeltaRDKB);
1119  }
1120  else if (fCandidateType == kD0toKpi) {
1121 
1122  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1123  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1124  fHistAlphaDpiR = new TH2F("fHistAlphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1125  fHistAlphaDKR = new TH2F("fHistAlphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1126 
1127  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1128  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1129 
1130  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1131  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1132  fHistDeltaRDpiR = new TH2F("fHistDeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1133  fHistDeltaRDKR = new TH2F("fHistDeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1134 
1135  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1136  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1137 
1138  fOutput->Add(fHistAlphaDpiS);
1139  fOutput->Add(fHistAlphaDKS);
1140  fOutput->Add(fHistAlphaDpiR);
1141  fOutput->Add(fHistAlphaDKR);
1142  fOutput->Add(fHistAlphaDpiB);
1143  fOutput->Add(fHistAlphaDKB);
1144 
1145  fOutput->Add(fHistDeltaRDpiS);
1146  fOutput->Add(fHistDeltaRDKS);
1147  fOutput->Add(fHistDeltaRDpiR);
1148  fOutput->Add(fHistDeltaRDKR);
1149  fOutput->Add(fHistDeltaRDpiB);
1150  fOutput->Add(fHistDeltaRDKB);
1151  }
1152  }
1153 
1154  return kTRUE;
1155 }
1156 
1157 //_______________________________________________________________________________
1158 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const
1159 {
1160  // Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1161 
1162  if (!p1 || !p2) return -1;
1163 
1164  Double_t phi1 = p1->Phi();
1165  Double_t eta1 = p1->Eta();
1166  Double_t phi2 = p2->Phi();
1167  Double_t eta2 = p2->Eta();
1168 
1169  Double_t dPhi = phi1 - phi2;
1170  if(dPhi < -(TMath::Pi())/2) dPhi = dPhi + TMath::TwoPi();
1171  if(dPhi > (3*(TMath::Pi()))/2) dPhi = dPhi - TMath::TwoPi();
1172 
1173  Double_t dEta = eta1 - eta2;
1174 
1175  Double_t deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1176 
1177  return deltaR;
1178 }
1179 
1180 //_______________________________________________________________________________
1182 {
1183  //
1184  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1185  //
1186 
1187  if (!tracks) return;
1188 
1189  TObjArray allDaughters(10);
1190  allDaughters.SetOwner(kFALSE);
1191 
1192  TIter next(coll);
1193  AliEmcalParticle* emcpart = 0;
1194  while ((emcpart = static_cast<AliEmcalParticle*>(next()))) {
1195  AliAODRecoDecay* reco = dynamic_cast<AliAODRecoDecay*>(emcpart->GetTrack());
1196  AliDebug(2, Form("Found a D meson candidtate with pT = %.3f, eta = %.3f, phi = %.3f", reco->Pt(), reco->Eta(), reco->Phi()));
1197  if (reco) AddDaughters(reco, allDaughters);
1198  }
1199 
1200  tracks->ResetCurrentID();
1201  AliVTrack* track = 0;
1202  Int_t n = coll->GetEntriesFast();
1203  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
1204  if (allDaughters.Remove(track) == 0) {
1205  new ((*coll)[n]) AliEmcalParticle(track);
1206  n++;
1207  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1208  }
1209  else {
1210  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1211  }
1212  }
1213 }
1214 
1215 //_______________________________________________________________________________
1217 {
1218  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1219 
1220  Int_t n = cand->GetNDaughters();
1221 
1222  //Printf("AddDaughters: the number of dauhters is %d", n);
1223 
1224  Int_t ntot = 0;
1225  Double_t pt = 0;
1226  for (Int_t i = 0; i < n; i++) {
1227  AliVTrack* track = dynamic_cast<AliVTrack*>(cand->GetDaughter(i));
1228  if (!track) continue;
1229 
1230  AliAODRecoDecay* cand2 = dynamic_cast<AliAODRecoDecay*>(track);
1231 
1232  if (cand2) {
1233  //Printf("Daughter pT = %.3f --> ", track->Pt());
1234  pt += AddDaughters(cand2, daughters);
1235  }
1236  else {
1237  if (!track->InheritsFrom("AliAODTrack")) {
1238  Printf("Warning: One of the daughters is not of type 'AliAODTrack' nor 'AliAODRecoDecay'.");
1239  continue;
1240  }
1241  //Printf("Daughter pT = %.3f", track->Pt());
1242  daughters.AddLast(track);
1243  pt += track->Pt();
1244  ntot++;
1245  }
1246  }
1247 
1248  //Printf("Total pt of the daughters = %.3f", pt);
1249 
1250  return pt;
1251 }
1252 
1253 //_________________________________________________________________________________________________
1254 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, TClonesArray* mcArray)
1255 {
1256  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1257 
1258  if (!mcArray) return -1;
1259 
1260  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1261  if (labDau0 < 0) return -1;
1262 
1263  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(mcArray->At(labDau0));
1264  return CheckOrigin(part, mcArray);
1265 }
1266 
1267 //_________________________________________________________________________________________________
1268 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, AliStack* stack)
1269 {
1270  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1271 
1272  if (!stack) return -1;
1273 
1274  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1275  if (labDau0 < 0) return -1;
1276 
1277  return CheckOrigin(labDau0, stack);
1278 }
1279 
1280 //_________________________________________________________________________________________________
1281 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODMCParticle* part, TClonesArray* mcArray)
1282 {
1283  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1284 
1285  if (!part) return -1;
1286  if (!mcArray) return -1;
1287 
1288  Int_t pdgGranma = 0;
1289  Int_t mother = part->GetMother();
1290  Int_t istep = 0;
1291  Int_t abspdgGranma = 0;
1292  Bool_t isFromB = kFALSE;
1293  Bool_t isQuarkFound = kFALSE;
1294 
1295  while (mother >= 0) {
1296  istep++;
1297  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1298  if (mcGranma >= 0) {
1299  pdgGranma = mcGranma->GetPdgCode();
1300  abspdgGranma = TMath::Abs(pdgGranma);
1301  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1302  isFromB = kTRUE;
1303  }
1304 
1305  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1306  mother = mcGranma->GetMother();
1307  }
1308  else {
1309  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1310  break;
1311  }
1312  }
1313 
1314  if (isQuarkFound) {
1315  if (isFromB) {
1316  return kFromBottom;
1317  }
1318  else {
1319  return kFromCharm;
1320  }
1321  }
1322  else {
1323  return kQuarkNotFound;
1324  }
1325 }
1326 
1327 //_________________________________________________________________________________________________
1329 {
1330  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1331 
1332  if (!stack) return -1;
1333 
1334  TParticle* part = stack->Particle(ipart);
1335  if (!part) return -1;
1336 
1337  Int_t pdgGranma = 0;
1338  Int_t mother = part->GetFirstMother();
1339  Int_t istep = 0;
1340  Int_t abspdgGranma = 0;
1341  Bool_t isFromB = kFALSE;
1342  Bool_t isQuarkFound = kFALSE;
1343 
1344  while (mother >= 0) {
1345  istep++;
1346  TParticle* mcGranma = stack->Particle(mother);
1347  if (mcGranma >= 0) {
1348  pdgGranma = mcGranma->GetPdgCode();
1349  abspdgGranma = TMath::Abs(pdgGranma);
1350  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1351  isFromB = kTRUE;
1352  }
1353 
1354  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1355  mother = mcGranma->GetFirstMother();
1356  }
1357  else {
1358  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1359  break;
1360  }
1361  }
1362 
1363  if (isQuarkFound) {
1364  if (isFromB) {
1365  return kFromBottom;
1366  }
1367  else {
1368  return kFromCharm;
1369  }
1370  }
1371  else {
1372  return kQuarkNotFound;
1373  }
1374 }
1375 
1376 //_________________________________________________________________________________________________
1377 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckDecayChannel(AliAODMCParticle* part, TClonesArray* mcArray)
1378 {
1379  // Determine the decay channel
1380 
1381  if (!part) return -1;
1382  if (!mcArray) return -1;
1383 
1385 
1386  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1387 
1388  if (part->GetNDaughters() == 2) {
1389 
1390  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1391  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1392 
1393  if (!d1 || !d2) {
1394  return decay;
1395  }
1396 
1397  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1398  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1399 
1400  if (absPdgPart == 421) { // D0 -> K pi
1401 
1402  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1403  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1404  decay = kDecayD0toKpi;
1405  }
1406  }
1407 
1408  if (absPdgPart == 413) { // D* -> D0 pi
1409 
1410  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1411  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1412  if (D0decay == kDecayD0toKpi) {
1413  decay = kDecayDStartoKpipi;
1414  }
1415  }
1416 
1417  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1418  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1419  if (D0decay == kDecayD0toKpi) {
1420  decay = kDecayDStartoKpipi;
1421  }
1422  }
1423  }
1424  }
1425 
1426  return decay;
1427 }
1428 
1429 //_________________________________________________________________________________________________
1431 {
1432  // Determine the decay channel
1433 
1434  if (!stack) return -1;
1435 
1436  TParticle* part = stack->Particle(ipart);
1437 
1438  if (!part) return -1;
1439 
1441 
1442  if (part->GetNDaughters() == 2) {
1443 
1444  Int_t id1 = part->GetDaughter(0);
1445  Int_t id2 = part->GetDaughter(1);
1446 
1447  TParticle* d1 = stack->Particle(id1);
1448  TParticle* d2 = stack->Particle(id2);
1449 
1450  if (!d1 || !d2) {
1451  return decay;
1452  }
1453 
1454  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1455  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1456 
1457 
1458  if (part->GetPdgCode() == 421) { // D0 -> K pi
1459 
1460  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1461  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1462  decay = kDecayD0toKpi;
1463  }
1464  }
1465 
1466  if (part->GetPdgCode() == 413) { // D* -> D0 pi
1467 
1468  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1469  Int_t D0decay = CheckDecayChannel(id1, stack);
1470  if (D0decay == kDecayD0toKpi) {
1471  decay = kDecayDStartoKpipi;
1472  }
1473  }
1474 
1475  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1476  Int_t D0decay = CheckDecayChannel(id2, stack);
1477  if (D0decay == kDecayD0toKpi) {
1478  decay = kDecayDStartoKpipi;
1479  }
1480  }
1481  }
1482  }
1483 
1484  return decay;
1485 
1486 }
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.
Double_t mass
void FillDStarMCTruthKinHistos(AliAODRecoCascadeHF *dstar, Int_t, Int_t isDstar)
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
Container for particles within the EMCAL framework.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Int_t fNCand
contains bkg candidates selected by AliRDHFCuts and the rest of the event tracks
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
const Double_t ptmin
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
ClassImp(AliAnalysisTaskSEDmesonsFilterCJ) AliAnalysisTaskSEDmesonsFilterCJ
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:276
Double_t InvMassD0() const
TH1 * fHistStat
number of selected side-band D candidates already added to fSideBandArray
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:234
Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:296
Int_t PtBin(Double_t pt) const
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t nptbins
TClonesArray * fCombinedDmesons
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...