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