AliPhysics  857879d (857879d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 fCandidateType(),
65 fCorrelationMethod(),
66 fPDGmother(),
67 fNProngs(),
68 fPDGdaughters(),
69 fBranchName(),
70 fCuts(0),
71 fMinMass(),
72 fMaxMass(),
73 fCandidateArray(0),
74 fSideBandArray(0),
75 fAnalyseDBkg(kFALSE),
76 fBuildRM(kFALSE),
77 fBuildRMEff(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 fCandidateType(),
109 fCorrelationMethod(),
110 fPDGmother(),
111 fNProngs(),
112 fPDGdaughters(),
113 fBranchName(),
114 fCuts(0),
115 fMinMass(),
116 fMaxMass(),
117 fCandidateArray(0),
118 fSideBandArray(0),
119 fAnalyseDBkg(kFALSE),
120 fBuildRM(kFALSE),
121 fBuildRMEff(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 fillRM[8] = {zRec,JetPtRec,Drec->Pt(),etaRec,zGen,MCjet->Pt(),Dgen->Pt(),etaGen};
570  fResponseMatrix->Fill(fillRM,1.);
571  }
572 }
573 
574 //_______________________________________________________________________________
575 
577 {
578 
579  if(!MCjet) AliDebug(2, "No Generated Level Jet Found!");
580 
581  Int_t pdgMeson = 413;
582  if (fCandidateType == kD0toKpi) pdgMeson = 421;
583 
584  AliVParticle *Dgen = MCjet->GetFlavourTrack();
585  Double_t zGen = Z(Dgen,MCjet,0);
586  Double_t JetPtGen = MCjet->Pt();
587  Double_t JetEtaGen = MCjet->Eta();
588  Double_t DPtGen = Dgen->Pt();
589  Double_t DEtaGen = Dgen->Eta();
590 
591  Double_t zRec = -999;
592  Double_t JetPtRec = -999;
593  Double_t DPtRec = -999;
594  Double_t JetEtaRec = -999;
595 
596  AliJetContainer* JetContRec = GetJetContainer(0);
597  if(JetContRec){
598  AliEmcalJet* jet;
599  GetHFJet(jet,kFALSE);
600 
601  if (jet){
602 
603  AliVParticle *Drec = jet->GetFlavourTrack();
604 
605  Double_t rho = 0;
606  if(!fUsePythia){
607  rho = JetContRec->GetRhoVal();
608  if(fLocalRho) rho = fLocalRho->GetLocalVal(jet->Phi(),JetContRec->GetJetRadius(),JetContRec->GetRhoVal());
609  }
610  zRec = 0;
611  if(rho>0) zRec = Z(Drec,jet,rho);
612  else zRec = Z(Drec,jet);
613  JetPtRec = jet->Pt() - rho*jet->Area();
614  JetEtaRec = jet->Eta();
615 
616  DPtRec = Drec->Pt();
617 
618  Bool_t bDInEMCalAcc=InEMCalAcceptance(Drec);
619  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
620 
622  {
623  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Drec;
624  FillHistogramsD0JetCorr(dzero,zRec,Drec->Pt(),JetPtRec,kFALSE,bDInEMCalAcc,bJetInEMCalAcc,aodEvent);
625  }
626  else if(fCandidateType==kDstartoKpipi)
627  {
629  FillHistogramsDstarJetCorr(dstar,zRec,Drec->Pt(),JetPtRec,kFALSE,bDInEMCalAcc,bJetInEMCalAcc);
630  }
631 
632  } // if HF reco jet
633  } // if jet cont reco
634 
635  Double_t fillRM[8] = {zRec,JetPtRec,DPtRec,JetEtaRec,zGen,JetPtGen,DPtGen,JetEtaGen};
636  fResponseMatrix->Fill(fillRM,1.);
637 
638 }
639 
640 //_______________________________________________________________________________
641 
643 {
644 
645  AliVParticle *Dmeson = jet->GetFlavourTrack(0);
646  Double_t JetPtCorr = jet->Pt() - rho*jet->Area();
647  Double_t z = 0;
648  if(rho>0) z = Z(Dmeson,jet,rho);
649  else z = Z(Dmeson,jet);
650 
651  if(IsBkg==kFALSE && fBuildRM==kTRUE) CreateResponseMatrix(jet);
652 
653  Bool_t bDInEMCalAcc=InEMCalAcceptance(Dmeson);
654  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
655 
657  {
658  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Dmeson;
659  FillHistogramsD0JetCorr(dzero,z,Dmeson->Pt(),JetPtCorr,IsBkg,bDInEMCalAcc,bJetInEMCalAcc,aodEvent);
660  }
661  else if(fCandidateType==kDstartoKpipi)
662  {
664  FillHistogramsDstarJetCorr(dstar,z,Dmeson->Pt(),JetPtCorr,IsBkg,bDInEMCalAcc,bJetInEMCalAcc);
665  }
666 
667 
668 }
669 
671 {
672  AliJetContainer* JetCont = 0;
673 
674  if(!IsBkg) JetCont = GetJetContainer(0);
675  else JetCont = GetJetContainer(1);
676 
677  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
678 
679  Bool_t JetIsHF = kFALSE;
680 
681  JetCont->ResetCurrentID();
682  while ((jet = JetCont->GetNextJet()))
683  {
684  UInt_t rejectionReason = 0;
685  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
686  if(!OKjet) continue;
687 
688  Int_t JetTag = AliEmcalJet::kD0;
689  TString recoDecayClassName("AliAODRecoDecayHF2Prong");
691  {
692  JetTag = AliEmcalJet::kDStar;
693  recoDecayClassName = "AliAODRecoCascadeHF";
694  }
695  //loop on jet particles
696  Int_t ntrjet= jet->GetNumberOfTracks();
697  for(Int_t itrk=0;itrk<ntrjet;itrk++)
698  {
699  AliVParticle* jetTrk=jet->TrackAt(itrk,ParticlesCont->GetArray());
700  if(!jetTrk) continue;
701  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
702  if(emcpart) jetTrk = emcpart->GetTrack();
703  if(strcmp(jetTrk->ClassName(),recoDecayClassName)==0)
704  {
705  JetIsHF = kTRUE;
706  if(!IsBkg) jet->AddFlavourTag(JetTag);
707  jet->AddFlavourTrack(jetTrk);
708  break;
709  }
710  } //end loop on jet tracks
711  if(JetIsHF) break;
712  } //end of jet loop
713 
714  if(!JetIsHF) jet = 0;
715 
716 }
717 //_______________________________________________________________________________
718 
720 {
721  Bool_t HFMCjet = kFALSE;
722 
723  AliJetContainer* mcjets = 0;
724 
725  if(!fAnalyseDBkg) mcjets = GetJetContainer(1);
726  else mcjets = GetJetContainer(2);
727 
728  AliParticleContainer *ParticlesCont = mcjets->GetParticleContainer();
729  mcjets->ResetCurrentID();
730 
731  Int_t njet=0;
732 
733  while ((mcjet = mcjets->GetNextAcceptJet()))
734  {
735  njet++;
736  //loop on jet particles
737  Int_t ntrjet= mcjet->GetNumberOfTracks();
738 
739  for(Int_t itrk=0;itrk<ntrjet;itrk++)
740  {
741  AliAODMCParticle* jetTrk=(AliAODMCParticle*)mcjet->TrackAt(itrk,ParticlesCont->GetArray());
742 
743  if(TMath::Abs(jetTrk->GetPdgCode())==fPDGmother)
744  {
745  HFMCjet=kTRUE;
746  mcjet->AddFlavourTrack(jetTrk);
747  break;
748  }
749  } //end loop on jet tracks
750  if(HFMCjet==kTRUE) break;
751  }
752 
753  if(!HFMCjet) mcjet=0;
754 }
755 
756 //_______________________________________________________________________________
757 
759 {
760  // The Terminate() function is the last function to be called during
761  // a query. It always runs on the client, it can be used to present
762  // the results graphically or save the results to file.
763 
764  Info("Terminate"," terminate");
765  AliAnalysisTaskSE::Terminate();
766 
767  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
768  if (!fOutput) {
769  printf("ERROR: fOutput not available\n");
770  return;
771  }
772 }
773 
774 //_______________________________________________________________________________
775 
777  Float_t mass=0;
778  Int_t abspdg=TMath::Abs(pdg);
779 
780  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
781  // compute the Delta mass for the D*
783  Float_t mass1=0;
784  mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
785  mass = mass-mass1;
786  }
787 
788  fMinMass = mass-range;
789  fMaxMass = mass+range;
790 
791  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
792  if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
793 }
794 
795 //_______________________________________________________________________________
796 
798  if(uplimit>lowlimit)
799  {
800  fMinMass = lowlimit;
801  fMaxMass = uplimit;
802  }
803  else{
804  printf("Error! Lower limit larger than upper limit!\n");
805  fMinMass = uplimit - uplimit*0.2;
806  fMaxMass = uplimit;
807  }
808 }
809 
810 //_______________________________________________________________________________
811 
813  if(nptbins>30) {
814  AliInfo("Maximum number of bins allowed is 30!");
815  return kFALSE;
816  }
817  if(!width) return kFALSE;
818  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
819  return kTRUE;
820 }
821 
822 //_______________________________________________________________________________
823 
825 
826  Double_t p[3],pj[3];
827  Bool_t okpp=part->PxPyPz(p);
828  Bool_t okpj=jet->PxPyPz(pj);
829 
830  if(!okpp || !okpj)
831  {
832  printf("Problems getting momenta\n");
833  return -999;
834  }
835 
836  //Background Subtraction
837  //It corrects the each component of the jet momentum for Z calculation
838 
839  pj[0] = jet->Px() - jet->Area()*(rho*TMath::Cos(jet->AreaPhi()));
840  pj[1] = jet->Py() - jet->Area()*(rho*TMath::Sin(jet->AreaPhi()));
841  pj[2] = jet->Pz() - jet->Area()*(rho*TMath::SinH(jet->AreaEta()));
842 
843  return Z(p,pj);
844 }
845 //_______________________________________________________________________________
846 
848 
849  Double_t p[3],pj[3];
850  Bool_t okpp=part->PxPyPz(p);
851  Bool_t okpj=jet->PxPyPz(pj);
852 
853  if(!okpp || !okpj)
854  {
855  printf("Problems getting momenta\n");
856  return -999;
857  }
858 
859  return Z(p,pj);
860 }
861 //_______________________________________________________________________________
863 
864  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
865  Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
866  return z;
867 }
868 
869 
870 //_______________________________________________________________________________
872 
873  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
874  Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
875  return z;
876 }
877 
878 //_______________________________________________________________________________
879 
881 
882  // Statistics
883  Int_t nbins=8;
884  if(fUseMCInfo) nbins+=2;
885 
886  fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
887  fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
888  fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
889  fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
890  fhstat->GetXaxis()->SetBinLabel(4,"N jets");
891  fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
892  fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
893  fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
894  fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
895  if(fUseMCInfo) {
896  fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
897  fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
898  fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
899  fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
900  }
901  fhstat->SetNdivisions(1);
902  fOutput->Add(fhstat);
903 
904  fhCentDjet=new TH1F("hCentDjet","Centrality",100,0,100);
905  fOutput->Add(fhCentDjet);
906 
907  const Int_t nbinsmass=300;
908  const Int_t nbinspttrack=500;
909  const Int_t nbinsptjet=200;
910  const Int_t nbinsptD=100;
911  const Int_t nbinsphi=200;
912  const Int_t nbinseta=100;
913 
914  //binning for THnSparse
915  const Int_t nbinsSpsmass=120;
916  const Int_t nbinsSpsptjet=200;
917  const Int_t nbinsSpsptD=100;
918  const Int_t nbinsSpsz=160;
919  const Int_t nbinsSpsy=150;
920 
921  const Float_t pttracklims[2]={0.,200.};
922  const Float_t ptjetlims[2]={-50.,150.};
923  const Float_t ptDlims[2]={0.,50.};
924  const Float_t zlims[2]={-1.2,2};
925  const Float_t philims[2]={0.,6.3};
926  const Float_t etalims[2]={-1.5,1.5};
927 
928  // jet related fistograms
929 
930  fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
931  fhPhiJetTrks->Sumw2();
932  fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
933  fhEtaJetTrks->Sumw2();
934  fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
935  fhPtJetTrks->Sumw2();
936 
937  fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
938  fhPhiJet->Sumw2();
939  fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
940  fhEtaJet->Sumw2();
941  fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
942  fhPtJet->Sumw2();
943 
944  fOutput->Add(fhPhiJetTrks);
945  fOutput->Add(fhEtaJetTrks);
946  fOutput->Add(fhPtJetTrks);
947  fOutput->Add(fhPhiJet);
948  fOutput->Add(fhEtaJet);
949  fOutput->Add(fhPtJet);
950 
952  {
953  if(fAnalyseDBkg){
954  fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
955  fhDiffSideBand->SetStats(kTRUE);
956  fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
957  fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
958  fhDiffSideBand->Sumw2();
959  fOutput->Add(fhDiffSideBand);
960  }
961 
962  fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
963  fhPtPion->SetStats(kTRUE);
964  fhPtPion->GetXaxis()->SetTitle("GeV/c");
965  fhPtPion->GetYaxis()->SetTitle("Entries");
966  fhPtPion->Sumw2();
967  fOutput->Add(fhPtPion);
968 
969  }
970 
971 
972  fhInvMassptD = new TH2F("hInvMassptD","D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
973  fhInvMassptD->SetStats(kTRUE);
974  fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
975  fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
976  fhInvMassptD->Sumw2();
977 
978  fOutput->Add(fhInvMassptD);
979 
980  if(fUseMCInfo){
981  fhInvMassptDbg = new TH2F("hInvMassptDbg","Background D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
982  fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
983  fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
984  fhInvMassptDbg->Sumw2();
985  fOutput->Add(fhInvMassptDbg);
986 
987  }
988 
989  fhsDphiz=0x0; //definition below according to the switches
990 
991  if(fUseMCInfo){
992  AliInfo("Creating a 8 axes container (MB background candidates)");
993  const Int_t nAxis=8;
994  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,nbinsSpsy, 2, 2, 2};
995  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,etalims[0], -0.5,-0.5,-0.5};
996  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,etalims[1], 1.5, 1.5 , 1.5};
997  fNAxesBigSparse=nAxis;
998  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., y^{D}, Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
999 
1000  } else{
1001  AliInfo("Creating a 5 axes container");
1002  const Int_t nAxis=5;
1003  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,nbinsSpsy};
1004  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,etalims[0]};
1005  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,etalims[1]};
1006  fNAxesBigSparse=nAxis;
1007 
1008  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., y^{D}", nAxis, nbinsSparse, minSparse, maxSparse);
1009  }
1010 
1011  if(!fhsDphiz) AliFatal("No THnSparse created");
1012  fhsDphiz->Sumw2();
1013 
1014  fOutput->Add(fhsDphiz);
1015 
1016  if(fBuildRM == kTRUE || fBuildRMEff == kTRUE)
1017  {
1018  const Int_t nRMAxis=8;
1019  const Int_t RMnbins[nRMAxis] = {nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsy,nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsy};
1020  const Double_t RMmin[nRMAxis]={zlims[0],ptjetlims[0],ptDlims[0],etalims[0],zlims[0],ptjetlims[0],ptDlims[0],etalims[0]};
1021  const Double_t RMmax[nRMAxis]={zlims[1],ptjetlims[1],ptDlims[1],etalims[1],zlims[1],ptjetlims[1],ptDlims[1],etalims[1]};
1022  fResponseMatrix = new THnSparseF("ResponseMatrix","z, jet pt, D pt, #eta^{jet}: Rec and Gen",nRMAxis,RMnbins,RMmin,RMmax);
1023  fResponseMatrix->Sumw2();
1024  fOutput->Add(fResponseMatrix);
1025  }
1026 
1027  PostData(1, fOutput);
1028 
1029  return kTRUE;
1030 }
1031 
1032 //_______________________________________________________________________________
1033 
1035 
1036 
1037  Double_t masses[2]={0.,0.};
1038  Int_t pdgdaughtersD0[2]={211,321};//pi,K
1039  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1040  Int_t pdgMeson = 413;
1041  if (fCandidateType == kD0toKpi) pdgMeson = 421;
1042 
1043  masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1044  masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1045 
1046  Double_t *point=0x0;
1047 
1048  if(!fAnalyseDBkg)
1049  {
1050  point=new Double_t[5];
1051  point[0]=z;
1052  point[1]=ptj;
1053  point[2]=ptD;
1054  point[3]=masses[0];
1055  point[4]=candidate->Y(pdgMeson);
1056 
1057  }
1058  else
1059  {
1060  point=new Double_t[8];
1061  point[0]=z;
1062  point[1]=ptj;
1063  point[2]=ptD;
1064  point[3]=masses[0];
1065  point[4]=candidate->Y(pdgMeson);
1066  point[5]=static_cast<Double_t>(IsBkg ? 1 : 0);
1067  point[6]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1068  point[7]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1069 
1070  }
1071  Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1072  if(isselected==1 || isselected==3)
1073  {
1074 
1075  fhInvMassptD->Fill(masses[0],ptD);
1076  fhsDphiz->Fill(point,1.);
1077 
1078  }
1079  if(isselected>=2)
1080  {
1081  fhInvMassptD->Fill(masses[1],ptD);
1082  point[3]=masses[1];
1083  fhsDphiz->Fill(point,1.);
1084 
1085  }
1086  delete[] point;
1087 }
1088 
1089 //_______________________________________________________________________________
1090 
1092 
1093  AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1094  Double_t deltamass= dstar->DeltaInvMass();
1095  Int_t pdgMeson = 413;
1096  if (fCandidateType == kD0toKpi) pdgMeson = 421;
1097 
1098  fhPtPion->Fill(softpi->Pt());
1099 
1100  Double_t *point=0x0;
1101  if(!fAnalyseDBkg)
1102  {
1103  point=new Double_t[5];
1104  point[0]=z;
1105  point[1]=ptj;
1106  point[2]=ptD;
1107  point[3]=deltamass;
1108  point[4]=dstar->Y(pdgMeson);
1109  }
1110  else
1111  {
1112  point=new Double_t[8];
1113  point[0]=z;
1114  point[1]=ptj;
1115  point[2]=ptD;
1116  point[3]=deltamass;
1117  point[4]=dstar->Y(pdgMeson);
1118  point[5]=static_cast<Double_t>(IsBkg ? 1 : 0);
1119  point[6]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1120  point[7]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1121  }
1122 
1123  if(!point){
1124  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1125  return;
1126  }
1127 
1128 
1129  fhInvMassptD->Fill(deltamass,ptD);
1130  fhsDphiz->Fill(point,1.);
1131  delete[] point;
1132 }
1133 
1134 //_______________________________________________________________________________
1135 
1137 
1138  Double_t pdgmass=0;
1139  if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1140  if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1141 
1142  Double_t *point=0x0;
1143 
1144  if(fNAxesBigSparse==5){
1145  point=new Double_t[5];
1146  point[0]=z;
1147  point[1]=ptjet;
1148  point[2]=ptD;
1149  point[3]=pdgmass;
1150  point[4]=yD;
1151  }
1152  if(fNAxesBigSparse==8){
1153  point=new Double_t[8];
1154  point[0]=z;
1155  point[1]=ptjet;
1156  point[2]=ptD;
1157  point[3]=pdgmass;
1158  point[4]=yD;
1159  point[5]=1;
1160  point[6]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1161  point[7]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1162  }
1163 
1164  if(!point){
1165  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1166  return;
1167  }
1168 
1169 
1170  point[3]=pdgmass;
1171  fhsDphiz->Fill(point,1.);
1172 
1173  delete[] point;
1174 }
1175 
1176 //_______________________________________________________________________________
1177 
1179  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1180  //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1181  if(!p1 || !p2) return -1;
1182 
1183  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1184 
1185  if (rho>0)
1186  {
1187  Double_t pj[3];
1188  Bool_t okpj=p1->PxPyPz(pj);
1189  if(!okpj){
1190  printf("Problems getting momenta\n");
1191  return -999;
1192  }
1193  pj[0] = p1->Px() - p1->Area()*(rho*TMath::Cos(p1->AreaPhi()));
1194  pj[1] = p1->Py() - p1->Area()*(rho*TMath::Sin(p1->AreaPhi()));
1195  pj[2] = p1->Pz() - p1->Area()*(rho*TMath::SinH(p1->AreaEta()));
1196  //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1197  //[0,pi] if py > 0 and
1198  //[pi,2pi] if py < 0
1199  phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1200  if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1201  eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1202  }
1203 
1204  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1205 
1206  Double_t dPhi=TMath::Abs(phi1-phi2);
1207  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1208 
1209 
1210  Double_t dEta=eta1-eta2;
1211  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1212  return deltaR;
1213 
1214 }
1215 //_______________________________________________________________________________
1217  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1218  if(!p1 || !p2) return -1;
1219 
1220  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1221 
1222  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1223 
1224  Double_t dPhi=TMath::Abs(phi1-phi2);
1225  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1226 
1227 
1228  Double_t dEta=eta1-eta2;
1229  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1230  return deltaR;
1231 
1232 }
1233 
1234 //_______________________________________________________________________________
1235 
1237 
1238  Double_t ptD=candDstar->Pt();
1239  Int_t bin = fCuts->PtBin(ptD);
1240  if (bin < 0){
1241  // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1242  bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1243  return -1; // out of bounds
1244  }
1245 
1246  Double_t invM=candDstar->InvMassD0();
1247  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1248 
1249  Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1250  Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1251 
1252  if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1253  else return 0;
1254 
1255 }
1256 
1257 //_______________________________________________________________________________
1258 
1260  //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1261 
1262  Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1263  Bool_t binEMCal=kTRUE;
1264  Double_t phi=vpart->Phi(), eta=vpart->Eta();
1265  if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1266  if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1267  return binEMCal;
1268 
1269 
1270 }
Int_t pdg
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:246
Double_t Area() const
Definition: AliEmcalJet.h:123
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
void FillHistogramsMCGenDJetCorr(Double_t z, Double_t ptD, Double_t ptjet, Double_t yD, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
Double_t DeltaInvMass() const
Jet is tagged to contain a D* meson.
Definition: AliEmcalJet.h:78
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Z(AliVParticle *part, AliEmcalJet *jet, Double_t rho) const
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
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()
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF *dstar, Double_t z, Double_t ptD, Double_t ptj, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
AliVTrack * GetTrack() const
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:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
AliParticleContainer * GetParticleContainer() const
Double_t AreaEta() const
Definition: AliEmcalJet.h:125
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
void FillHistogramsD0JetCorr(AliAODRecoDecayHF *candidate, Double_t z, Double_t ptD, Double_t ptj, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent *aodEvent)
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()
ClassImp(AliAnalysisTaskFlavourJetCorrelations) AliAnalysisTaskFlavourJetCorrelations
Double_t AreaPhi() const
Definition: AliEmcalJet.h:126
Double_t Pt() const
Definition: AliEmcalJet.h:102
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:79
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:281
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Bool_t PxPyPz(Double_t p[3]) const
Definition: AliEmcalJet.h:104
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
Double_t Pz() const
Definition: AliEmcalJet.h:101
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:239
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...