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