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