AliPhysics  e6c8d43 (e6c8d43)
AliAnalysisTaskFlavourJetCorrelations.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 // Analysis Taks for reconstructed particle correlation
18 // (first implementation done for D mesons) with jets
19 // (use the so called Emcal framework)
20 //
21 //-----------------------------------------------------------------------
22 // Authors:
23 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
24 // S. Antônio (University of São Paulo) antonio.silva@cern.ch
25 // A. Grelli (Utrecht University) a.grelli@uu.nl
26 // X. Zhang (LBNL) XMZhang@lbl.gov
27 // B. Trzeciak (Utrecht University) barbara.antonina.trzeciak@cern.ch
28 //-----------------------------------------------------------------------
29 
30 #include <TDatabasePDG.h>
31 #include <TParticle.h>
32 #include "TROOT.h"
33 #include <THnSparse.h>
34 #include <TSystem.h>
35 #include <TObjectTable.h>
36 #include "AliMultSelection.h"
37 
39 #include "AliAODHandler.h"
40 #include "AliAnalysisManager.h"
41 #include "AliEmcalJet.h"
42 #include "AliJetContainer.h"
43 #include "AliAODRecoDecay.h"
44 #include "AliAODRecoCascadeHF.h"
46 #include "AliAODMCParticle.h"
47 #include "AliRDHFCutsD0toKpi.h"
49 #include "AliParticleContainer.h"
50 #include "AliEmcalParticle.h"
51 #include "AliLocalRhoParameter.h"
53 
55 
56 
57 //_______________________________________________________________________________
58 
61 fUseMCInfo(kTRUE),
62 fUseReco(kTRUE),
63 fUsePythia(kFALSE),
64 fBuildRM(kFALSE),
65 fBuildRMEff(kFALSE),
66 fCandidateType(),
67 fCorrelationMethod(),
68 fPDGmother(),
69 fNProngs(),
70 fPDGdaughters(),
71 fBranchName(),
72 fCuts(0),
73 fMinMass(),
74 fMaxMass(),
75 fCandidateArray(0),
76 fSideBandArray(0),
77 fAnalyseDBkg(kFALSE),
78 fNAxesBigSparse(9),
79 fUseCandArray(kFALSE),
80 fUseSBArray(kFALSE),
81 fhstat(),
82 fhCentDjet(),
83 fhPtJetTrks(),
84 fhPhiJetTrks(),
85 fhEtaJetTrks(),
86 fhPtJet(),
87 fhPhiJet(),
88 fhEtaJet(),
89 fhInvMassptD(),
90 fhDiffSideBand(),
91 fhInvMassptDbg(),
92 fhPtPion(),
93 fhsDphiz(),
94 fResponseMatrix()
95 
96 {
97  //
98  // Default ctor
99 }
100 
101 //_______________________________________________________________________________
102 
104 AliAnalysisTaskEmcalJet(name,kTRUE),
105 fUseMCInfo(kTRUE),
106 fUseReco(kTRUE),
107 fUsePythia(kFALSE),
108 fBuildRM(kFALSE),
109 fBuildRMEff(kFALSE),
110 fCandidateType(),
111 fCorrelationMethod(),
112 fPDGmother(),
113 fNProngs(),
114 fPDGdaughters(),
115 fBranchName(),
116 fCuts(0),
117 fMinMass(),
118 fMaxMass(),
119 fCandidateArray(0),
120 fSideBandArray(0),
121 fAnalyseDBkg(kFALSE),
122 fNAxesBigSparse(9),
123 fUseCandArray(kFALSE),
124 fUseSBArray(kFALSE),
125 fhstat(),
126 fhCentDjet(),
127 fhPtJetTrks(),
128 fhPhiJetTrks(),
129 fhEtaJetTrks(),
130 fhPtJet(),
131 fhPhiJet(),
132 fhEtaJet(),
133 fhInvMassptD(),
134 fhDiffSideBand(),
135 fhInvMassptDbg(),
136 fhPtPion(),
137 fhsDphiz(),
138 fResponseMatrix()
139 {
140  //
141  // Constructor. Initialization of Inputs and Outputs
142  //
143 
144  Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
145  fCuts=cuts;
146  fCandidateType=candtype;
147  const Int_t nptbins=fCuts->GetNPtBins();
148  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};
149 
150  switch(fCandidateType){
151  case 0:
152  fPDGmother=421;
153  fNProngs=2;
154  fPDGdaughters[0]=211;//pi
155  fPDGdaughters[1]=321;//K
156  fPDGdaughters[2]=0; //empty
157  fPDGdaughters[3]=0; //empty
158  fBranchName="D0toKpi";
159  break;
160  case 1:
161  fPDGmother=413;
162  fNProngs=3;
163  fPDGdaughters[1]=211;//pi soft
164  fPDGdaughters[0]=421;//D0
165  fPDGdaughters[2]=211;//pi fromD0
166  fPDGdaughters[3]=321; // K from D0
167  fBranchName="Dstar";
168 
169  if(nptbins<20){
170  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
171  } else {
172  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
173  }
174  break;
175  default:
176  printf("%d not accepted!!\n",fCandidateType);
177  break;
178  }
179 
182  if(fUseCandArray) DefineInput(1, TClonesArray::Class());
183  if(fUseSBArray) DefineInput(2, TClonesArray::Class());
184 
185  DefineOutput(1,TList::Class()); // histos
186  DefineOutput(2,AliRDHFCuts::Class()); // my cuts
187 
188 }
189 
190 //_______________________________________________________________________________
191 
193  //
194  // destructor
195  //
196 
197  Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
198 
199  delete fCuts;
200 
201 }
202 
203 //_______________________________________________________________________________
204 
206  //
207  // Initialization
208  //
209 
210  if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
211 
212  switch(fCandidateType){
213  case 0:
214  {
215  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
216  copyfCuts->SetName("AnalysisCutsDzero");
217  // Post the data
218  PostData(2,copyfCuts);
219  }
220  break;
221  case 1:
222  {
223  AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
224  copyfCuts->SetName("AnalysisCutsDStar");
225  // Post the cuts
226  PostData(2,copyfCuts);
227  }
228  break;
229  default:
230  return;
231  }
232 
233  return;
234 }
235 
236 //_______________________________________________________________________________
237 
239  // output
240  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
242 
243  // define histograms
244  // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
246  PostData(1,fOutput);
247 
248  return;
249 }
250 
251 //_______________________________________________________________________________
252 
254 {
255  // user exec from AliAnalysisTaskEmcal is used
256 
257  // Load the event
258  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
259 
260  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
261  if (matchingAODdeltaAODlevel<=0) {
262  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
263  return kFALSE;
264  }
265 
266  TClonesArray *arrayDStartoD0pi=0;
267 
268  if (!aodEvent && AODEvent() && IsStandardAOD()) {
269 
270  // In case there is an AOD handler writing a standard AOD, use the AOD
271  // event in memory rather than the input (ESD) event.
272  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
273 
274  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
275  // have to taken from the AOD event hold by the AliAODExtension
276  AliAODHandler* aodHandler = (AliAODHandler*)
277  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
278  if(aodHandler->GetExtensions()) {
279  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
280  AliAODEvent *aodFromExt = ext->GetAOD();
281  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
282  }
283  } else if(aodEvent){
284  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
285  }
286 
287  if (!arrayDStartoD0pi) {
288  AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
289  // return;
290  } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
291 
292  TClonesArray* mcArray = 0x0;
293  if (fUseMCInfo) { //not used at the moment,uncomment return if you use
294  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
295  if (!mcArray) {
296  printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
297  }
298  }
299 
300  //D meson candidates. Also background if is MC
301 
302  if(fUseCandArray)
303  {
304  fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
305  if (!fCandidateArray) return kFALSE;
306  for(Int_t icand=0; icand<fCandidateArray->GetEntriesFast(); icand++)
307  {
308  fhstat->Fill(2);
309  }
310  }
311  if(fUseSBArray)
312  {
313  fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
314  if (!fSideBandArray) return kFALSE;
315  }
316 
317  fhstat->Fill(0);
318 
319  // fix for temporary bug in ESDfilter
320  // the AODs with null vertex pointer didn't pass the PhysSel
321  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
322 
323  //Event selection
324  Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
325  TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
326  if(!iseventselected) return kFALSE;
327 
328  fhstat->Fill(1);
329  if(aodEvent->GetRunNumber()>200000)
330  {
331  Float_t lPercentile = 300;
332  AliMultSelection *MultSelection = 0x0;
333  MultSelection = (AliMultSelection*)aodEvent->FindListObject("MultSelection");
334  if(!MultSelection) {
335  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
336  AliWarning("AliMultSelection object not found!");
337  }else{
338  lPercentile = MultSelection->GetMultiplicityPercentile("V0M");
339  }
340 
341  fhCentDjet->Fill(lPercentile);
342  }
343  else fhCentDjet->Fill(fCent);
344 
345 
346 // for MC response matrix of efficiency studies, fMultCand option only
347 if(fUseMCInfo && fBuildRMEff){
348 
349  AliJetContainer* mcjets = 0;
350  if(!fAnalyseDBkg) mcjets = GetJetContainer(1);
351  else mcjets = GetJetContainer(2);
352  if(!mcjets) return kFALSE;
353 
354  AliParticleContainer *MCParticlesCont = mcjets->GetParticleContainer();
355 
356  mcjets->ResetCurrentID();
357  AliEmcalJet* jet=0;
358 
359  while ((jet = mcjets->GetNextJet()))
360  {
361  UInt_t rejectionReason = 0;
362  Bool_t OKjet = mcjets->AcceptJet(jet, rejectionReason);
363  if(!OKjet) {
364  fhstat->Fill(5);
365  continue;
366  }
367 
368  fhstat->Fill(3); //Jet accepted
369  fhPhiJet->Fill(jet->Phi());
370  fhEtaJet->Fill(jet->Eta());
371  fhPtJet->Fill(jet->Pt());
372 
373  Int_t ntrjet= jet->GetNumberOfTracks();
374 
375  for(Int_t itrk=0;itrk<ntrjet;itrk++)
376  {
377  AliAODMCParticle* jetTrk=(AliAODMCParticle*)jet->TrackAt(itrk,MCParticlesCont->GetArray());
378  if (!jetTrk) continue;
379  fhPtJetTrks->Fill(jetTrk->Pt());
380  fhPhiJetTrks->Fill(jetTrk->Phi());
381  fhEtaJetTrks->Fill(jetTrk->Eta());
382  } //end loop on jet tracks
383 
384  } // end loop on mc jets
385 
386  // Get HF accepted MC jet
387  AliEmcalJet* MCjet = 0;
388  FindMCJet(MCjet);
389  if(!MCjet) return kFALSE;
390  //if( TMath::Abs(MCjet->Eta()) > (0.9 - mcjets->GetJetRadius()) ) return kFALSE;
391 
393  {
394  if(fBuildRMEff==kTRUE) CreateMCResponseMatrix(MCjet, aodEvent);
395  }
396  /* the other method not enabled for now
397  * else if(fCorrelationMethod==kAngular)
398  {
399  if(fCandidateArray->GetEntriesFast()>0) AngularCorrelationMethod(kFALSE,aodEvent);
400  if(fAnalyseDBkg==kTRUE && fSideBandArray->GetEntriesFast()>0) AngularCorrelationMethod(kTRUE,aodEvent);
401  }*/
402 
403 }
404 else {
405 
406  AliJetContainer* JetCont = GetJetContainer(0);
407  if(!JetCont) return kFALSE;
408  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
409 
410  AliJetContainer* JetContSB = 0;
411  AliParticleContainer *ParticlesContSB = 0;
412  if(fAnalyseDBkg)
413  {
414  JetContSB = GetJetContainer(1);
415  if(!JetContSB) return kFALSE;
416  ParticlesContSB = JetContSB->GetParticleContainer();
417  }
418 
419  //Distribution of all particles in the event
420  Int_t ntrarr=ParticlesCont->GetNParticles();
421  for(Int_t i=0;i<ntrarr;i++)
422  {
423  AliVParticle* jetTrk= ParticlesCont->GetParticle(i);
424  if (!jetTrk) continue;
425  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
426  if (emcpart) jetTrk = emcpart->GetTrack();
427  fhPtJetTrks->Fill(jetTrk->Pt());
428  fhPhiJetTrks->Fill(jetTrk->Phi());
429  fhEtaJetTrks->Fill(jetTrk->Eta());
430  }
431 
432  JetCont->ResetCurrentID();
433  AliEmcalJet* jet=0;
434 
435  while ((jet = JetCont->GetNextJet()))
436  {
437  UInt_t rejectionReason = 0;
438  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
439  if(!OKjet) {
440  fhstat->Fill(5);
441  continue;
442  }
443 
444  Double_t JetPtCorr = 0;
445  if(fUseMCInfo && fUsePythia){
446  JetPtCorr = jet->Pt();
447  }
448  else {
449  JetPtCorr = jet->Pt() - jet->Area()*JetCont->GetRhoVal(); //background subtraction
450  if(fLocalRho)
451  {
452  JetPtCorr = jet->Pt() - jet->Area()*fLocalRho->GetLocalVal(jet->Phi(),JetCont->GetJetRadius(),JetCont->GetRhoVal()); //
453  }
454  }
455  fhstat->Fill(3); //Jet accepted
456  fhPhiJet->Fill(jet->Phi());
457  fhEtaJet->Fill(jet->Eta());
458  fhPtJet->Fill(JetPtCorr);
459  }
460 
461  if(ParticlesCont->GetNParticles()>0) fhstat->Fill(2);
462 
464  {
465  if(ParticlesCont->GetNParticles()>0) ConstituentCorrelationMethod(kFALSE,aodEvent);
466  if(fAnalyseDBkg==kTRUE && ParticlesContSB->GetNParticles()>0) ConstituentCorrelationMethod(kTRUE,aodEvent);
467  }
468 
469  else if(fCorrelationMethod==kAngular)
470  {
471  if(fCandidateArray->GetEntriesFast()>0) AngularCorrelationMethod(kFALSE,aodEvent);
472  if(fAnalyseDBkg==kTRUE && fSideBandArray->GetEntriesFast()>0) AngularCorrelationMethod(kTRUE,aodEvent);
473  }
474 }
475 
476 
477  PostData(1,fOutput);
478  return kTRUE;
479 }
481 {
482 
483  AliJetContainer* JetCont = 0;
484 
485  if(!IsBkg) JetCont = GetJetContainer(0);
486  else JetCont = GetJetContainer(1);
487 
488  Double_t rho = 0;
489  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
490 
491  AliEmcalJet *jet;
492 
493  GetHFJet(jet,IsBkg);
494 
495  if(jet)
496  {
497  if(fLocalRho) rho = fLocalRho->GetLocalVal(jet->Phi(),JetCont->GetJetRadius(),JetCont->GetRhoVal());
498  FillDJetHistograms(jet,rho,IsBkg,aodEvent);
499  }
500 
501 }
503 {
504  AliJetContainer* JetCont = GetJetContainer(0);
505 
506  Int_t ncand = 0;
507  if(!IsBkg) ncand = fCandidateArray->GetEntriesFast();
508  else ncand = fSideBandArray->GetEntriesFast();
509 
510  Double_t rho = 0;
511  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
512 
513  for(Int_t icand = 0; icand<ncand; icand++)
514  {
515  AliVParticle* charm=0x0;
516  if(!IsBkg) charm = (AliVParticle*)fCandidateArray->At(icand);
517  else charm = (AliVParticle*)fSideBandArray->At(icand);
518  if(!charm) continue;
519 
520  Int_t JetTag = AliEmcalJet::kD0;
522  //loop over jets
523  JetCont->ResetCurrentID();
524  AliEmcalJet* jet=0;
525  while ((jet = JetCont->GetNextJet()))
526  {
527  UInt_t rejectionReason = 0;
528  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
529  if(!OKjet) continue;
530 
531  if(fLocalRho) rho = fLocalRho->GetLocalVal(jet->Phi(),JetCont->GetJetRadius(),JetCont->GetRhoVal());
532 
533  if(DeltaR(jet,charm,rho)<JetCont->GetJetRadius() && CheckDeltaR(jet,charm)<JetCont->GetJetRadius())
534  {
535  if(!IsBkg) jet->AddFlavourTag(JetTag);
536  jet->AddFlavourTrack(charm);
537  FillDJetHistograms(jet,rho,IsBkg,aodEvent);
538  }
539  }
540  }
541 
542 }
543 
544 //_______________________________________________________________________________
545 
547 {
548  AliJetContainer* JetContRec = GetJetContainer(0);
549 
550  Double_t rho = 0;
551  if(!fUsePythia) rho = JetContRec->GetRhoVal();
552 
553  AliEmcalJet* MCjet = 0;
554 
555  FindMCJet(MCjet);
556 
557  if(!jet) AliDebug(2, "No Reconstructed Level Jet Found!");
558  else if(!MCjet) AliDebug(2, "No Generated Level Jet Found!");
559  else
560  {
561  if(fLocalRho) rho = fLocalRho->GetLocalVal(jet->Phi(),JetContRec->GetJetRadius(),JetContRec->GetRhoVal());
562  AliVParticle *Drec = jet->GetFlavourTrack();
563  AliVParticle *Dgen = MCjet->GetFlavourTrack();
564  Double_t zRec = Z(Drec,jet,rho);
565  Double_t zGen = Z(Dgen,MCjet,0);
566  Double_t JetPtRec = jet->Pt() - rho*jet->Area();
567  Double_t etaRec = jet->Eta();
568  Double_t etaGen = MCjet->Eta();
569  Double_t DYRec = -999;
570  AliAODMCParticle* DTrk=(AliAODMCParticle*)Dgen;
571  Double_t DYGen = DTrk->Y();
572 
573  Int_t pdgMeson = 413;
574  if (fCandidateType == kD0toKpi) pdgMeson = 421;
575 
577  {
578  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Drec;
579  DYRec = dzero->Y(pdgMeson);
580 
581  }
582  else if(fCandidateType==kDstartoKpipi)
583  {
585  DYRec = dstar->Y(pdgMeson);
586 
587  }
588 
589  Double_t fillRM[10] = {zRec,JetPtRec,Drec->Pt(),DYRec,etaRec,zGen,MCjet->Pt(),Dgen->Pt(),DYGen,etaGen};
590  fResponseMatrix->Fill(fillRM,1.);
591  }
592 }
593 
594 //_______________________________________________________________________________
595 
597 {
598  if(!MCjet) AliDebug(2, "No Generated Level Jet Found!");
599 
600  AliVParticle *Dgen = MCjet->GetFlavourTrack();
601  Double_t zGen = Z(Dgen,MCjet,0);
602  Double_t JetPtGen = MCjet->Pt();
603  Double_t JetEtaGen = MCjet->Eta();
604  Double_t DPtGen = Dgen->Pt();
605  AliAODMCParticle* DTrk = (AliAODMCParticle*)Dgen;
606  Double_t DYGen = DTrk->Y();
607  Int_t pdg = DTrk->GetPdgCode();
608 
609  Double_t zRec = -999;
610  Double_t JetPtRec = -999;
611  Double_t JetEtaRec = -999;
612  Double_t DPtRec = -999;
613  Double_t DYRec = -999;
614 
615  AliJetContainer* JetContRec = GetJetContainer(0);
616  if(JetContRec){
617  AliEmcalJet* jet;
618  GetHFJet(jet,kFALSE);
619 
620  if (jet){
621  AliVParticle *Drec = jet->GetFlavourTrack();
622 
623  Double_t rho = 0;
624  if(!fUsePythia){
625  rho = JetContRec->GetRhoVal();
626  if(fLocalRho) rho = fLocalRho->GetLocalVal(jet->Phi(),JetContRec->GetJetRadius(),JetContRec->GetRhoVal());
627  }
628  zRec = 0;
629  if(rho>0) zRec = Z(Drec,jet,rho);
630  else zRec = Z(Drec,jet);
631  JetPtRec = jet->Pt() - rho*jet->Area();
632  JetEtaRec = jet->Eta();
633  DPtRec = Drec->Pt();
634 
635  Int_t pdgMeson = 413;
636  if (fCandidateType == kD0toKpi) pdgMeson = 421;
637 
639  {
640  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Drec;
641  DYRec = dzero->Y(pdgMeson);
642 
643  }
644  else if(fCandidateType==kDstartoKpipi)
645  {
647  DYRec = dstar->Y(pdgMeson);
648 
649  }
650 
651  Bool_t bDInEMCalAcc=InEMCalAcceptance(Drec);
652  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
653 
655  {
656  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Drec;
657  FillHistogramsD0JetCorr(dzero,zRec,Drec->Pt(),JetPtRec,JetEtaRec,kFALSE,bDInEMCalAcc,bJetInEMCalAcc,aodEvent,pdg);
658  }
659  else if(fCandidateType==kDstartoKpipi)
660  {
662  FillHistogramsDstarJetCorr(dstar,zRec,Drec->Pt(),JetPtRec,JetEtaRec,kFALSE,bDInEMCalAcc,bJetInEMCalAcc);
663  }
664 
665  } // if HF reco jet
666  } // if jet cont reco
667 
668  Double_t fillRM[10] = {zRec,JetPtRec,DPtRec,DYRec,JetEtaRec,zGen,JetPtGen,DPtGen,DYGen,JetEtaGen};
669  fResponseMatrix->Fill(fillRM,1.);
670 
671 }
672 
673 //_______________________________________________________________________________
674 
676 {
677 
678  AliVParticle *Dmeson = jet->GetFlavourTrack(0);
679  Double_t JetPtCorr = jet->Pt() - rho*jet->Area();
680  Double_t JetEtaRec = jet->Eta();
681  Double_t z = 0;
682  if(rho>0) z = Z(Dmeson,jet,rho);
683  else z = Z(Dmeson,jet);
684 
685  if(IsBkg==kFALSE && fBuildRM==kTRUE) CreateResponseMatrix(jet);
686 
687  Bool_t bDInEMCalAcc=InEMCalAcceptance(Dmeson);
688  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
689 
691  {
692  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Dmeson;
693  FillHistogramsD0JetCorr(dzero,z,Dmeson->Pt(),JetPtCorr,JetEtaRec,IsBkg,bDInEMCalAcc,bJetInEMCalAcc,aodEvent,-999);
694  }
695  else if(fCandidateType==kDstartoKpipi)
696  {
698  FillHistogramsDstarJetCorr(dstar,z,Dmeson->Pt(),JetPtCorr,JetEtaRec,IsBkg,bDInEMCalAcc,bJetInEMCalAcc);
699  }
700 
701 
702 }
703 
705 {
706  AliJetContainer* JetCont = 0;
707 
708  if(!IsBkg) JetCont = GetJetContainer(0);
709  else JetCont = GetJetContainer(1);
710 
711  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
712 
713  Bool_t JetIsHF = kFALSE;
714 
715  JetCont->ResetCurrentID();
716  while ((jet = JetCont->GetNextJet()))
717  {
718  UInt_t rejectionReason = 0;
719  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
720  if(!OKjet) continue;
721 
722  Int_t JetTag = AliEmcalJet::kD0;
723  TString recoDecayClassName("AliAODRecoDecayHF2Prong");
725  {
726  JetTag = AliEmcalJet::kDStar;
727  recoDecayClassName = "AliAODRecoCascadeHF";
728  }
729  //loop on jet particles
730  Int_t ntrjet= jet->GetNumberOfTracks();
731  for(Int_t itrk=0;itrk<ntrjet;itrk++)
732  {
733  AliVParticle* jetTrk=jet->TrackAt(itrk,ParticlesCont->GetArray());
734  if(!jetTrk) continue;
735  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
736  if(emcpart) jetTrk = emcpart->GetTrack();
737  if(strcmp(jetTrk->ClassName(),recoDecayClassName)==0)
738  {
739  JetIsHF = kTRUE;
740  if(!IsBkg) jet->AddFlavourTag(JetTag);
741  jet->AddFlavourTrack(jetTrk);
742  break;
743  }
744  } //end loop on jet tracks
745  if(JetIsHF) break;
746  } //end of jet loop
747 
748  if(!JetIsHF) jet = 0;
749 
750 }
751 //_______________________________________________________________________________
752 
754 {
755  Bool_t HFMCjet = kFALSE;
756 
757  AliJetContainer* mcjets = 0;
758 
759  if(!fAnalyseDBkg) mcjets = GetJetContainer(1);
760  else mcjets = GetJetContainer(2);
761 
762  AliParticleContainer *ParticlesCont = mcjets->GetParticleContainer();
763  mcjets->ResetCurrentID();
764 
765  Int_t njet=0;
766 
767  while ((mcjet = mcjets->GetNextAcceptJet()))
768  {
769  njet++;
770  //loop on jet particles
771  Int_t ntrjet= mcjet->GetNumberOfTracks();
772 
773  for(Int_t itrk=0;itrk<ntrjet;itrk++)
774  {
775  AliAODMCParticle* jetTrk=(AliAODMCParticle*)mcjet->TrackAt(itrk,ParticlesCont->GetArray());
776 
777  if(TMath::Abs(jetTrk->GetPdgCode())==fPDGmother)
778  {
779  HFMCjet=kTRUE;
780  mcjet->AddFlavourTrack(jetTrk);
781  break;
782  }
783  } //end loop on jet tracks
784  if(HFMCjet==kTRUE) break;
785  }
786 
787  if(!HFMCjet) mcjet=0;
788 }
789 
790 //_______________________________________________________________________________
791 
793 {
794  // The Terminate() function is the last function to be called during
795  // a query. It always runs on the client, it can be used to present
796  // the results graphically or save the results to file.
797 
798  Info("Terminate"," terminate");
799  AliAnalysisTaskSE::Terminate();
800 
801  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
802  if (!fOutput) {
803  printf("ERROR: fOutput not available\n");
804  return;
805  }
806 }
807 
808 //_______________________________________________________________________________
809 
811  Float_t mass=0;
812  Int_t abspdg=TMath::Abs(pdg);
813 
814  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
815  // compute the Delta mass for the D*
817  Float_t mass1=0;
818  mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
819  mass = mass-mass1;
820  }
821 
822  fMinMass = mass-range;
823  fMaxMass = mass+range;
824 
825  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
826  if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
827 }
828 
829 //_______________________________________________________________________________
830 
832  if(uplimit>lowlimit)
833  {
834  fMinMass = lowlimit;
835  fMaxMass = uplimit;
836  }
837  else{
838  printf("Error! Lower limit larger than upper limit!\n");
839  fMinMass = uplimit - uplimit*0.2;
840  fMaxMass = uplimit;
841  }
842 }
843 
844 //_______________________________________________________________________________
845 
847  if(nptbins>30) {
848  AliInfo("Maximum number of bins allowed is 30!");
849  return kFALSE;
850  }
851  if(!width) return kFALSE;
852  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
853  return kTRUE;
854 }
855 
856 //_______________________________________________________________________________
857 
859 
860  Double_t p[3],pj[3];
861  Bool_t okpp=part->PxPyPz(p);
862  Bool_t okpj=jet->PxPyPz(pj);
863 
864  if(!okpp || !okpj)
865  {
866  printf("Problems getting momenta\n");
867  return -999;
868  }
869 
870  //Background Subtraction
871  //It corrects the each component of the jet momentum for Z calculation
872 
873  pj[0] = jet->Px() - jet->Area()*(rho*TMath::Cos(jet->AreaPhi()));
874  pj[1] = jet->Py() - jet->Area()*(rho*TMath::Sin(jet->AreaPhi()));
875  pj[2] = jet->Pz() - jet->Area()*(rho*TMath::SinH(jet->AreaEta()));
876 
877  return Z(p,pj);
878 }
879 //_______________________________________________________________________________
880 
882 
883  Double_t p[3],pj[3];
884  Bool_t okpp=part->PxPyPz(p);
885  Bool_t okpj=jet->PxPyPz(pj);
886 
887  if(!okpp || !okpj)
888  {
889  printf("Problems getting momenta\n");
890  return -999;
891  }
892 
893  return Z(p,pj);
894 }
895 //_______________________________________________________________________________
897 
898  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
899  Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
900  return z;
901 }
902 
903 
904 //_______________________________________________________________________________
906 
907  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
908  Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
909  return z;
910 }
911 
912 //_______________________________________________________________________________
913 
915 
916  // Statistics
917  Int_t nbins=8;
918  if(fUseMCInfo) nbins+=2;
919 
920  fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
921  fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
922  fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
923  fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
924  fhstat->GetXaxis()->SetBinLabel(4,"N jets");
925  fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
926  fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
927  fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
928  fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
929  if(fUseMCInfo) {
930  fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
931  fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
932  fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
933  fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
934  }
935  fhstat->SetNdivisions(1);
936  fOutput->Add(fhstat);
937 
938  fhCentDjet=new TH1F("hCentDjet","Centrality",100,0,100);
939  fOutput->Add(fhCentDjet);
940 
941  const Int_t nbinsmass=300;
942  const Int_t nbinspttrack=500;
943  const Int_t nbinsptjet=200;
944  const Int_t nbinsptD=100;
945  const Int_t nbinsphi=200;
946  const Int_t nbinseta=100;
947 
948  //binning for THnSparse
949  const Int_t nbinsSpsmass=120;
950  const Int_t nbinsSpsptjet=200;
951  const Int_t nbinsSpsptD=100;
952  const Int_t nbinsSpsz=160;
953  const Int_t nbinsSpsy=150;
954 
955  const Float_t pttracklims[2]={0.,200.};
956  const Float_t ptjetlims[2]={-50.,150.};
957  const Float_t ptDlims[2]={0.,50.};
958  const Float_t zlims[2]={-1.2,2};
959  const Float_t philims[2]={0.,6.3};
960  const Float_t etalims[2]={-1.5,1.5};
961 
962  // jet related fistograms
963 
964  fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
965  fhPhiJetTrks->Sumw2();
966  fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
967  fhEtaJetTrks->Sumw2();
968  fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
969  fhPtJetTrks->Sumw2();
970 
971  fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
972  fhPhiJet->Sumw2();
973  fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
974  fhEtaJet->Sumw2();
975  fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
976  fhPtJet->Sumw2();
977 
978  fOutput->Add(fhPhiJetTrks);
979  fOutput->Add(fhEtaJetTrks);
980  fOutput->Add(fhPtJetTrks);
981  fOutput->Add(fhPhiJet);
982  fOutput->Add(fhEtaJet);
983  fOutput->Add(fhPtJet);
984 
986  {
987  if(fAnalyseDBkg){
988  fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
989  fhDiffSideBand->SetStats(kTRUE);
990  fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
991  fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
992  fhDiffSideBand->Sumw2();
993  fOutput->Add(fhDiffSideBand);
994  }
995 
996  fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
997  fhPtPion->SetStats(kTRUE);
998  fhPtPion->GetXaxis()->SetTitle("GeV/c");
999  fhPtPion->GetYaxis()->SetTitle("Entries");
1000  fhPtPion->Sumw2();
1001  fOutput->Add(fhPtPion);
1002 
1003  }
1004 
1005 
1006  fhInvMassptD = new TH2F("hInvMassptD","D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
1007  fhInvMassptD->SetStats(kTRUE);
1008  fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1009  fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1010  fhInvMassptD->Sumw2();
1011 
1012  fOutput->Add(fhInvMassptD);
1013 
1014  if(fUseMCInfo){
1015  fhInvMassptDbg = new TH2F("hInvMassptDbg","Background D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
1016  fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1017  fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1018  fhInvMassptDbg->Sumw2();
1019  fOutput->Add(fhInvMassptDbg);
1020 
1021  }
1022 
1023  fhsDphiz=0x0; //definition below according to the switches
1024 
1025  if(fUseMCInfo){
1027  AliInfo("Creating a 11 axes container (MB background candidates)");
1028  const Int_t nAxis=11;
1029  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,nbinsSpsy,nbinsSpsy,nbinsSpsmass,nbinsSpsmass, 2, 2, 2};
1030  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,etalims[0],etalims[0],fMinMass,fMinMass, -0.5,-0.5,-0.5};
1031  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,etalims[1],etalims[1],fMaxMass,fMaxMass, 1.5, 1.5 , 1.5};
1032  fNAxesBigSparse=nAxis;
1033  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., y^{D}, #eta^{jet}, mass_{true}, mass_{refl}, Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1034 
1035  }
1036  else{
1037  AliInfo("Creating a 9 axes container (MB background candidates)");
1038  const Int_t nAxis=9;
1039  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,nbinsSpsy,nbinsSpsy, 2, 2, 2};
1040  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,etalims[0],etalims[0], -0.5,-0.5,-0.5};
1041  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,etalims[1],etalims[1], 1.5, 1.5 , 1.5};
1042  fNAxesBigSparse=nAxis;
1043  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., y^{D}, #eta^{jet}, Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1044  }
1045 
1046  } else{
1047  AliInfo("Creating a 6 axes container");
1048  const Int_t nAxis=6;
1049  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,nbinsSpsy,nbinsSpsy};
1050  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,etalims[0],etalims[0]};
1051  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,etalims[1],etalims[1]};
1052  fNAxesBigSparse=nAxis;
1053  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., y^{D}, #eta^{jet}", nAxis, nbinsSparse, minSparse, maxSparse);
1054  }
1055 
1056  if(!fhsDphiz) AliFatal("No THnSparse created");
1057  fhsDphiz->Sumw2();
1058 
1059  fOutput->Add(fhsDphiz);
1060 
1061  if(fBuildRM == kTRUE || fBuildRMEff == kTRUE)
1062  {
1063  const Int_t nRMAxis=10;
1064  const Int_t RMnbins[nRMAxis] = {nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsy,nbinsSpsy,nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsy,nbinsSpsy};
1065  const Double_t RMmin[nRMAxis]={zlims[0],ptjetlims[0],ptDlims[0],etalims[0],etalims[0],zlims[0],ptjetlims[0],ptDlims[0],etalims[0],etalims[0]};
1066  const Double_t RMmax[nRMAxis]={zlims[1],ptjetlims[1],ptDlims[1],etalims[1],etalims[1],zlims[1],ptjetlims[1],ptDlims[1],etalims[1],etalims[1]};
1067  fResponseMatrix = new THnSparseF("ResponseMatrix","z, p_{T}^{jet}, p_{T}^{D}, y^{D}, #eta^{jet}: Rec and Gen",nRMAxis,RMnbins,RMmin,RMmax);
1068  fResponseMatrix->Sumw2();
1069  fOutput->Add(fResponseMatrix);
1070  }
1071 
1072  PostData(1, fOutput);
1073 
1074  return kTRUE;
1075 }
1076 
1077 //_______________________________________________________________________________
1078 
1079 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t jetEta, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent, Int_t pdgTrue){
1080 
1081 
1082  Double_t masses[2]={0.,0.};
1083  Int_t pdgdaughtersD0[2]={211,321};//pi,K
1084  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1085  Int_t pdgMeson = 413;
1086  if (fCandidateType == kD0toKpi) pdgMeson = 421;
1087 
1088  masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1089  masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1090 
1091  Double_t *point=0x0;
1092 
1093  if(!fUseMCInfo)
1094  {
1095  point=new Double_t[6];
1096  point[0]=z;
1097  point[1]=ptj;
1098  point[2]=ptD;
1099  point[3]=masses[0];
1100  point[4]=candidate->Y(pdgMeson);
1101  point[5]=jetEta;
1102  }
1103  else
1104  {
1105  point=new Double_t[11];
1106  point[0]=z;
1107  point[1]=ptj;
1108  point[2]=ptD;
1109  point[3]=masses[0];
1110  point[4]=candidate->Y(pdgMeson);
1111  point[5]=jetEta;
1112  point[6]=-999;
1113  point[7]=-999;
1114  point[8]=static_cast<Double_t>(IsBkg ? 1 : 0);
1115  point[9]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1116  point[10]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1117 
1118  }
1119  Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1120 
1121  if(!fUseMCInfo) {
1122  if(isselected==1 || isselected==3)
1123  {
1124  fhInvMassptD->Fill(masses[0],ptD);
1125  point[3]=masses[0];
1126  fhsDphiz->Fill(point,1.);
1127  }
1128  if(isselected>=2)
1129  {
1130  fhInvMassptD->Fill(masses[1],ptD);
1131  point[3]=masses[1];
1132  fhsDphiz->Fill(point,1.);
1133  }
1134  }
1135  else {
1136  if(isselected==1 || isselected==3)
1137  {
1138  fhInvMassptD->Fill(masses[0],ptD);
1139  point[3]=masses[0];
1140 
1141  if(isselected==1) {
1142  point[6]=masses[0];
1143  }
1144  else if(isselected==3 && pdgTrue==421){
1145  point[6]=masses[0];
1146  point[7]=masses[1];
1147  }
1148  else if(isselected==3 && pdgTrue==-421){
1149  point[6]=masses[1];
1150  point[7]=masses[0];
1151  }
1152  fhsDphiz->Fill(point,1.);
1153  }
1154  if(isselected>=2)
1155  {
1156  fhInvMassptD->Fill(masses[1],ptD);
1157  point[3]=masses[1];
1158  if(isselected==2) {
1159  point[6]=masses[1];
1160  }
1161  fhsDphiz->Fill(point,1.);
1162  }
1163 
1164  }
1165 
1166  delete[] point;
1167 }
1168 
1169 //_______________________________________________________________________________
1170 
1172 
1173  AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1174  Double_t deltamass= dstar->DeltaInvMass();
1175  Int_t pdgMeson = 413;
1176  if (fCandidateType == kD0toKpi) pdgMeson = 421;
1177 
1178  fhPtPion->Fill(softpi->Pt());
1179 
1180  Double_t *point=0x0;
1181  if(!fAnalyseDBkg)
1182  {
1183  point=new Double_t[6];
1184  point[0]=z;
1185  point[1]=ptj;
1186  point[2]=ptD;
1187  point[3]=deltamass;
1188  point[4]=dstar->Y(pdgMeson);
1189  point[5]=jetEta;
1190  }
1191  else
1192  {
1193  point=new Double_t[9];
1194  point[0]=z;
1195  point[1]=ptj;
1196  point[2]=ptD;
1197  point[3]=deltamass;
1198  point[4]=dstar->Y(pdgMeson);
1199  point[5]=jetEta;
1200  point[6]=static_cast<Double_t>(IsBkg ? 1 : 0);
1201  point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1202  point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1203  }
1204 
1205  if(!point){
1206  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1207  return;
1208  }
1209 
1210 
1211  fhInvMassptD->Fill(deltamass,ptD);
1212  fhsDphiz->Fill(point,1.);
1213  delete[] point;
1214 }
1215 
1216 //_______________________________________________________________________________
1217 
1219 
1220  Double_t pdgmass=0;
1221  if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1222  if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1223 
1224  Double_t *point=0x0;
1225 
1226  if(fNAxesBigSparse==6){
1227  point=new Double_t[6];
1228  point[0]=z;
1229  point[1]=ptjet;
1230  point[2]=ptD;
1231  point[3]=pdgmass;
1232  point[4]=yD;
1233  point[5]=jetEta;
1234  }
1235  if(fNAxesBigSparse==9){
1236  point=new Double_t[9];
1237  point[0]=z;
1238  point[1]=ptjet;
1239  point[2]=ptD;
1240  point[3]=pdgmass;
1241  point[4]=yD;
1242  point[5]=jetEta;
1243  point[6]=1;
1244  point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1245  point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1246  }
1247 
1248  if(!point){
1249  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1250  return;
1251  }
1252 
1253 
1254  point[3]=pdgmass;
1255  fhsDphiz->Fill(point,1.);
1256 
1257  delete[] point;
1258 }
1259 
1260 //_______________________________________________________________________________
1261 
1263  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1264  //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1265  if(!p1 || !p2) return -1;
1266 
1267  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1268 
1269  if (rho>0)
1270  {
1271  Double_t pj[3];
1272  Bool_t okpj=p1->PxPyPz(pj);
1273  if(!okpj){
1274  printf("Problems getting momenta\n");
1275  return -999;
1276  }
1277  pj[0] = p1->Px() - p1->Area()*(rho*TMath::Cos(p1->AreaPhi()));
1278  pj[1] = p1->Py() - p1->Area()*(rho*TMath::Sin(p1->AreaPhi()));
1279  pj[2] = p1->Pz() - p1->Area()*(rho*TMath::SinH(p1->AreaEta()));
1280  //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1281  //[0,pi] if py > 0 and
1282  //[pi,2pi] if py < 0
1283  phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1284  if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1285  eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1286  }
1287 
1288  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1289 
1290  Double_t dPhi=TMath::Abs(phi1-phi2);
1291  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1292 
1293 
1294  Double_t dEta=eta1-eta2;
1295  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1296  return deltaR;
1297 
1298 }
1299 //_______________________________________________________________________________
1301  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1302  if(!p1 || !p2) return -1;
1303 
1304  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1305 
1306  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1307 
1308  Double_t dPhi=TMath::Abs(phi1-phi2);
1309  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1310 
1311 
1312  Double_t dEta=eta1-eta2;
1313  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1314  return deltaR;
1315 
1316 }
1317 
1318 //_______________________________________________________________________________
1319 
1321 
1322  Double_t ptD=candDstar->Pt();
1323  Int_t bin = fCuts->PtBin(ptD);
1324  if (bin < 0){
1325  // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1326  bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1327  return -1; // out of bounds
1328  }
1329 
1330  Double_t invM=candDstar->InvMassD0();
1331  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1332 
1333  Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1334  Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1335 
1336  if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1337  else return 0;
1338 
1339 }
1340 
1341 //_______________________________________________________________________________
1342 
1344  //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1345 
1346  Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1347  Bool_t binEMCal=kTRUE;
1348  Double_t phi=vpart->Phi(), eta=vpart->Eta();
1349  if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1350  if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1351  return binEMCal;
1352 
1353 
1354 }
Int_t pdg
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:351
Double_t Area() const
Definition: AliEmcalJet.h:130
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
Double_t DeltaInvMass() const
Jet is tagged to contain a D* meson.
Definition: AliEmcalJet.h:85
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Z(AliVParticle *part, AliEmcalJet *jet, Double_t rho) const
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Double_t mass
Int_t GetNParticles() const
Double_t GetLocalVal(Double_t phi, Double_t r, Double_t n) const
char Char_t
Definition: External.C:18
static Int_t CheckMatchingAODdeltaAODevents()
void FillHistogramsMCGenDJetCorr(Double_t z, Double_t ptD, Double_t ptjet, Double_t yD, Double_t jetEta, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
AliVTrack * GetTrack() const
void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF *dstar, Double_t z, Double_t ptD, Double_t ptj, Double_t jetEta, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
void AngularCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
AliVParticle * GetFlavourTrack(Int_t i=0) const
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
void FillHistogramsD0JetCorr(AliAODRecoDecayHF *candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t jetEta, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent *aodEvent, Int_t pdg)
AliParticleContainer * GetParticleContainer() const
Double_t AreaEta() const
Definition: AliEmcalJet.h:132
int Int_t
Definition: External.C:63
void FillDJetHistograms(AliEmcalJet *jet, Double_t rho, Bool_t IsBkg, AliAODEvent *aodEvent)
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void AddFlavourTrack(AliVParticle *hftrack)
virtual AliVParticle * GetParticle(Int_t i=-1) const
virtual 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.
AliAODTrack * GetBachelor() const
Double_t fCent
!event centrality
AliLocalRhoParameter * fLocalRho
! local event rho
Float_t CheckDeltaR(AliEmcalJet *p1, AliVParticle *p2) const
AliEmcalJet * GetNextAcceptJet()
Double_t AreaPhi() const
Definition: AliEmcalJet.h:133
Double_t Pt() const
Definition: AliEmcalJet.h:109
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
AliEmcalJet * GetNextJet()
Bool_t IsEventSelected(AliVEvent *event)
Jet is tagged to contain a D0 meson.
Definition: AliEmcalJet.h:86
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
void ConstituentCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
void CreateMCResponseMatrix(AliEmcalJet *MCjet, AliAODEvent *aodEvent)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Bool_t PxPyPz(Double_t p[3]) const
Definition: AliEmcalJet.h:111
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
Double_t Pz() const
Definition: AliEmcalJet.h:108
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:248
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
bool Bool_t
Definition: External.C:53
Int_t PtBin(Double_t pt) const
Float_t DeltaR(AliEmcalJet *p1, AliVParticle *p2, Double_t rho) const
Int_t nptbins
Container for jet within the EMCAL jet framework.
TClonesArray * fSideBandArray
contains candidates selected by AliRDHFCuts
Bool_t fAnalyseDBkg
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...