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