AliPhysics  a9863a5 (a9863a5)
 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 //-----------------------------------------------------------------------
28 
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
31 #include "TROOT.h"
32 #include <THnSparse.h>
33 #include <TSystem.h>
34 #include <TObjectTable.h>
35 
37 #include "AliAODHandler.h"
38 #include "AliAnalysisManager.h"
39 #include "AliEmcalJet.h"
40 #include "AliJetContainer.h"
41 #include "AliAODRecoDecay.h"
42 #include "AliAODRecoCascadeHF.h"
44 #include "AliAODMCParticle.h"
45 #include "AliRDHFCutsD0toKpi.h"
47 #include "AliParticleContainer.h"
48 #include "AliEmcalParticle.h"
49 
51 
52 
53 //_______________________________________________________________________________
54 
57 fUseMCInfo(kTRUE),
58 fUseReco(kTRUE),
59 fCandidateType(),
60 fCorrelationMethod(),
61 fPDGmother(),
62 fNProngs(),
63 fPDGdaughters(),
64 fBranchName(),
65 fCuts(0),
66 fMinMass(),
67 fMaxMass(),
68 fCandidateArray(0),
69 fSideBandArray(0),
70 fAnalyseDBkg(kFALSE),
71 fNAxesBigSparse(9),
72 fUseCandArray(kFALSE),
73 fUseSBArray(kFALSE),
74 fhstat(),
75 fhCentDjet(),
76 fhPtJetTrks(),
77 fhPhiJetTrks(),
78 fhEtaJetTrks(),
79 fhPtJet(),
80 fhPhiJet(),
81 fhEtaJet(),
82 fhInvMassptD(),
83 fhDiffSideBand(),
84 fhInvMassptDbg(),
85 fhPtPion(),
86 fhsDphiz(),
87 fResponseMatrix()
88 
89 {
90  //
91  // Default ctor
92 }
93 
94 //_______________________________________________________________________________
95 
97 AliAnalysisTaskEmcalJet(name,kTRUE),
98 fUseMCInfo(kTRUE),
99 fUseReco(kTRUE),
100 fCandidateType(),
101 fCorrelationMethod(),
102 fPDGmother(),
103 fNProngs(),
104 fPDGdaughters(),
105 fBranchName(),
106 fCuts(0),
107 fMinMass(),
108 fMaxMass(),
109 fCandidateArray(0),
110 fSideBandArray(0),
111 fAnalyseDBkg(kFALSE),
112 fNAxesBigSparse(9),
113 fUseCandArray(kFALSE),
114 fUseSBArray(kFALSE),
115 fhstat(),
116 fhCentDjet(),
117 fhPtJetTrks(),
118 fhPhiJetTrks(),
119 fhEtaJetTrks(),
120 fhPtJet(),
121 fhPhiJet(),
122 fhEtaJet(),
123 fhInvMassptD(),
124 fhDiffSideBand(),
125 fhInvMassptDbg(),
126 fhPtPion(),
127 fhsDphiz(),
128 fResponseMatrix()
129 {
130  //
131  // Constructor. Initialization of Inputs and Outputs
132  //
133 
134  Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
135  fCuts=cuts;
136  fCandidateType=candtype;
137  const Int_t nptbins=fCuts->GetNPtBins();
138  Float_t defaultSigmaD013[13]={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};
139 
140  switch(fCandidateType){
141  case 0:
142  fPDGmother=421;
143  fNProngs=2;
144  fPDGdaughters[0]=211;//pi
145  fPDGdaughters[1]=321;//K
146  fPDGdaughters[2]=0; //empty
147  fPDGdaughters[3]=0; //empty
148  fBranchName="D0toKpi";
149  break;
150  case 1:
151  fPDGmother=413;
152  fNProngs=3;
153  fPDGdaughters[1]=211;//pi soft
154  fPDGdaughters[0]=421;//D0
155  fPDGdaughters[2]=211;//pi fromD0
156  fPDGdaughters[3]=321; // K from D0
157  fBranchName="Dstar";
158 
159  if(nptbins<=13){
160  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
161  } else {
162  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
163  }
164  break;
165  default:
166  printf("%d not accepted!!\n",fCandidateType);
167  break;
168  }
169 
172  if(fUseCandArray) DefineInput(1, TClonesArray::Class());
173  if(fUseSBArray) DefineInput(2, TClonesArray::Class());
174 
175  DefineOutput(1,TList::Class()); // histos
176  DefineOutput(2,AliRDHFCuts::Class()); // my cuts
177 
178 }
179 
180 //_______________________________________________________________________________
181 
183  //
184  // destructor
185  //
186 
187  Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
188 
189  delete fCuts;
190 
191 }
192 
193 //_______________________________________________________________________________
194 
196  //
197  // Initialization
198  //
199 
200  if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
201 
202  switch(fCandidateType){
203  case 0:
204  {
205  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
206  copyfCuts->SetName("AnalysisCutsDzero");
207  // Post the data
208  PostData(2,copyfCuts);
209  }
210  break;
211  case 1:
212  {
213  AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
214  copyfCuts->SetName("AnalysisCutsDStar");
215  // Post the cuts
216  PostData(2,copyfCuts);
217  }
218  break;
219  default:
220  return;
221  }
222 
223  return;
224 }
225 
226 //_______________________________________________________________________________
227 
229  // output
230  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
232 
233  // define histograms
234  // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
236  PostData(1,fOutput);
237 
238  return;
239 }
240 
241 //_______________________________________________________________________________
242 
244 {
245  // user exec from AliAnalysisTaskEmcal is used
246 
247  // Load the event
248  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
249 
250  TClonesArray *arrayDStartoD0pi=0;
251 
252  if (!aodEvent && AODEvent() && IsStandardAOD()) {
253 
254  // In case there is an AOD handler writing a standard AOD, use the AOD
255  // event in memory rather than the input (ESD) event.
256  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
257 
258  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
259  // have to taken from the AOD event hold by the AliAODExtension
260  AliAODHandler* aodHandler = (AliAODHandler*)
261  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
262  if(aodHandler->GetExtensions()) {
263  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
264  AliAODEvent *aodFromExt = ext->GetAOD();
265  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
266  }
267  } else if(aodEvent){
268  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
269  }
270 
271  if (!arrayDStartoD0pi) {
272  AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
273  // return;
274  } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
275 
276  TClonesArray* mcArray = 0x0;
277  if (fUseMCInfo) { //not used at the moment,uncomment return if you use
278  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
279  if (!mcArray) {
280  printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
281  //return kFALSE;
282  }
283  }
284 
285  //D meson candidates. Also background if is MC
286 
287  if(fUseCandArray)
288  {
289  fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
290  if (!fCandidateArray) return kFALSE;
291  for(Int_t icand=0; icand<fCandidateArray->GetEntriesFast(); icand++)
292  {
293  fhstat->Fill(2);
294  }
295  }
296  if(fUseSBArray)
297  {
298  fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
299  if (!fSideBandArray) return kFALSE;
300  }
301 
302  fhstat->Fill(0);
303 
304  // fix for temporary bug in ESDfilter
305  // the AODs with null vertex pointer didn't pass the PhysSel
306  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
307 
308  //Event selection
309  Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
310  TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
311  if(!iseventselected) return kFALSE;
312 
313  fhstat->Fill(1);
314  fhCentDjet->Fill(fCent);
315 
316  AliJetContainer* JetCont = GetJetContainer(0);
317  if(!JetCont) return kFALSE;
318  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
319 
320  AliJetContainer* JetContSB = 0;
321  AliParticleContainer *ParticlesContSB = 0;
322  if(fAnalyseDBkg)
323  {
324  JetContSB = GetJetContainer(1);
325  if(!JetContSB) return kFALSE;
326  ParticlesContSB = JetContSB->GetParticleContainer();
327  }
328 
329  //Distribution of all particles in the event
330  Int_t ntrarr=ParticlesCont->GetNParticles();
331  for(Int_t i=0;i<ntrarr;i++)
332  {
333  AliVParticle* jetTrk= ParticlesCont->GetParticle(i);
334  if (!jetTrk) continue;
335  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
336  if (emcpart) jetTrk = emcpart->GetTrack();
337  fhPtJetTrks->Fill(jetTrk->Pt());
338  fhPhiJetTrks->Fill(jetTrk->Phi());
339  fhEtaJetTrks->Fill(jetTrk->Eta());
340  }
341 
342  JetCont->ResetCurrentID();
343  AliEmcalJet* jet=0;
344 
345  while ((jet = JetCont->GetNextJet()))
346  {
347  UInt_t rejectionReason = 0;
348  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
349  if(!OKjet) {
350  fhstat->Fill(5);
351  continue;
352  }
353 
354 
355  Double_t JetPtCorr = jet->Pt() - jet->Area()*JetCont->GetRhoVal(); //background subtraction
356  fhstat->Fill(3); //Jet accepted
357  fhPhiJet->Fill(jet->Phi());
358  fhEtaJet->Fill(jet->Eta());
359  fhPtJet->Fill(JetPtCorr);
360  }
361 
362  if(ParticlesCont->GetNParticles()>0) fhstat->Fill(2);
363 
365  {
366  if(ParticlesCont->GetNParticles()>0) ConstituentCorrelationMethod(kFALSE,aodEvent);
367  if(fAnalyseDBkg==kTRUE && ParticlesContSB->GetNParticles()>0) ConstituentCorrelationMethod(kTRUE,aodEvent);
368  }
369 
370  else if(fCorrelationMethod==kAngular)
371  {
372  if(fCandidateArray->GetEntriesFast()>0) AngularCorrelationMethod(kFALSE,aodEvent);
373  if(fAnalyseDBkg==kTRUE && fSideBandArray->GetEntriesFast()>0) AngularCorrelationMethod(kTRUE,aodEvent);
374  }
376  {
378  }
379 
380  PostData(1,fOutput);
381  return kTRUE;
382 }
384 {
385 
386  AliJetContainer* JetCont = 0;
387 
388  if(!IsBkg) JetCont = GetJetContainer(0);
389  else JetCont = GetJetContainer(1);
390 
391  Double_t rho = 0;
392  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
393 
394  AliEmcalJet *jet;
395 
396  GetHFJet(jet,IsBkg);
397 
398  if(jet) FillDJetHistograms(jet,rho,IsBkg,aodEvent);
399 
400 }
402 {
403  AliJetContainer* JetCont = GetJetContainer(0);
404 
405  Int_t ncand = 0;
406  if(!IsBkg) ncand = fCandidateArray->GetEntriesFast();
407  else ncand = fSideBandArray->GetEntriesFast();
408 
409  Double_t rho = 0;
410  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
411 
412  for(Int_t icand = 0; icand<ncand; icand++)
413  {
414  AliVParticle* charm=0x0;
415  if(!IsBkg) charm = (AliVParticle*)fCandidateArray->At(icand);
416  else charm = (AliVParticle*)fSideBandArray->At(icand);
417  if(!charm) continue;
418 
419  Int_t JetTag = AliEmcalJet::kD0;
421  //loop over jets
422  JetCont->ResetCurrentID();
423  AliEmcalJet* jet=0;
424  while ((jet = JetCont->GetNextJet()))
425  {
426  UInt_t rejectionReason = 0;
427  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
428  if(!OKjet) continue;
429 
430  if(DeltaR(jet,charm,rho)<JetCont->GetJetRadius() && CheckDeltaR(jet,charm)<JetCont->GetJetRadius())
431  {
432  if(!IsBkg) jet->AddFlavourTag(JetTag);
433  jet->AddFlavourTrack(charm);
434  FillDJetHistograms(jet,rho,IsBkg,aodEvent);
435  }
436  }
437  }
438 
439 }
440 
442 {
443  AliJetContainer* JetContRec = GetJetContainer(0);
444  //AliJetContainer* JetContGen = GetJetContainer(1);
445 
446  Double_t rho = JetContRec->GetRhoVal();
447 
448  AliEmcalJet *jetRec = 0;
449  AliEmcalJet *jetGen = 0;
450 
451  GetHFJet(jetRec,kFALSE);
452  GetHFJet(jetGen,kTRUE);
453  //GetHFJetParticleLevel(jetGen,kTRUE);
454 
455  if(!jetRec) AliDebug(2, "No Reconstructed Level Jet Found!");
456  else if(!jetGen) AliDebug(2, "No Generated Level Jet Found!");
457  else
458  {
459  AliVParticle *Drec = jetRec->GetFlavourTrack();
460  AliVParticle *Dgen = jetGen->GetFlavourTrack();
461  Double_t zRec = Z(Drec,jetRec,rho);
462  Double_t zGen = Z(Dgen,jetGen,0);
463  Double_t JetPtRec = jetRec->Pt() - rho*jetRec->Area();
464  Double_t fillRM[6] = {zRec,JetPtRec,Drec->Pt(),zGen,jetGen->Pt(),Dgen->Pt()};
465  fResponseMatrix->Fill(fillRM,1.);
466  }
467 }
468 
469 //_______________________________________________________________________________
470 
471 void AliAnalysisTaskFlavourJetCorrelations::FillDJetHistograms(AliEmcalJet* jet, Double_t rho, Bool_t IsBkg, AliAODEvent* aodEvent)
472 {
473 
474  AliVParticle *Dmeson = jet->GetFlavourTrack(0);
475  Double_t JetPtCorr = jet->Pt() - rho*jet->Area();
476  Double_t z = 0;
477  if(rho>0) z = Z(Dmeson,jet,rho);
478  else z = Z(Dmeson,jet);
479 
480  Bool_t bDInEMCalAcc=InEMCalAcceptance(Dmeson);
481  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
482 
484  {
485  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Dmeson;
486  FillHistogramsD0JetCorr(dzero,z,Dmeson->Pt(),JetPtCorr,jet->Phi(),IsBkg,bDInEMCalAcc,bJetInEMCalAcc,aodEvent);
487  }
488  else if(fCandidateType==kDstartoKpipi)
489  {
491  FillHistogramsDstarJetCorr(dstar,z,Dmeson->Pt(),JetPtCorr,jet->Phi(),IsBkg,bDInEMCalAcc,bJetInEMCalAcc);
492  }
493 
494 
495 }
496 
498 {
499  AliJetContainer* JetCont = 0;
500 
501  if(!IsBkg) JetCont = GetJetContainer(0);
502  else JetCont = GetJetContainer(1);
503 
504  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
505 
506  Bool_t JetIsHF = kFALSE;
507 
508  JetCont->ResetCurrentID();
509  while ((jet = JetCont->GetNextJet()))
510  {
511  UInt_t rejectionReason = 0;
512  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
513  if(!OKjet) continue;
514 
515  Int_t JetTag = AliEmcalJet::kD0;
516  TString recoDecayClassName("AliAODRecoDecayHF2Prong");
518  {
519  JetTag = AliEmcalJet::kDStar;
520  recoDecayClassName = "AliAODRecoCascadeHF";
521  }
522  //loop on jet particles
523  Int_t ntrjet= jet->GetNumberOfTracks();
524  for(Int_t itrk=0;itrk<ntrjet;itrk++)
525  {
526  AliVParticle* jetTrk=jet->TrackAt(itrk,ParticlesCont->GetArray());
527  if(!jetTrk) continue;
528  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
529  if(emcpart) jetTrk = emcpart->GetTrack();
530  if(strcmp(jetTrk->ClassName(),recoDecayClassName)==0)
531  {
532  JetIsHF = kTRUE;
533  if(!IsBkg) jet->AddFlavourTag(JetTag);
534  jet->AddFlavourTrack(jetTrk);
535  break;
536  }
537  } //end loop on jet tracks
538  if(JetIsHF) break;
539  } //end of jet loop
540 
541  if(!JetIsHF) jet = 0;
542 
543 }
544 
545 //_______________________________________________________________________________
546 
548 {
549  // The Terminate() function is the last function to be called during
550  // a query. It always runs on the client, it can be used to present
551  // the results graphically or save the results to file.
552 
553  Info("Terminate"," terminate");
554  AliAnalysisTaskSE::Terminate();
555 
556  fOutput = dynamic_cast<TList*> (GetOutputData(1));
557  if (!fOutput) {
558  printf("ERROR: fOutput not available\n");
559  return;
560  }
561 }
562 
563 //_______________________________________________________________________________
564 
566  Float_t mass=0;
567  Int_t abspdg=TMath::Abs(pdg);
568 
569  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
570  // compute the Delta mass for the D*
572  Float_t mass1=0;
573  mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
574  mass = mass-mass1;
575  }
576 
577  fMinMass = mass-range;
578  fMaxMass = mass+range;
579 
580  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
581  if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
582 }
583 
584 //_______________________________________________________________________________
585 
586 void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
587  if(uplimit>lowlimit)
588  {
589  fMinMass = lowlimit;
590  fMaxMass = uplimit;
591  }
592  else{
593  printf("Error! Lower limit larger than upper limit!\n");
594  fMinMass = uplimit - uplimit*0.2;
595  fMaxMass = uplimit;
596  }
597 }
598 
599 //_______________________________________________________________________________
600 
602  if(nptbins>30) {
603  AliInfo("Maximum number of bins allowed is 30!");
604  return kFALSE;
605  }
606  if(!width) return kFALSE;
607  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
608  return kTRUE;
609 }
610 
611 //_______________________________________________________________________________
612 
613 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Double_t rho) const{
614 
615  Double_t p[3],pj[3];
616  Bool_t okpp=part->PxPyPz(p);
617  Bool_t okpj=jet->PxPyPz(pj);
618 
619  if(!okpp || !okpj)
620  {
621  printf("Problems getting momenta\n");
622  return -999;
623  }
624 
625  //Background Subtraction
626  //It corrects the each component of the jet momentum for Z calculation
627 
628  pj[0] = jet->Px() - jet->Area()*(rho*TMath::Cos(jet->AreaPhi()));
629  pj[1] = jet->Py() - jet->Area()*(rho*TMath::Sin(jet->AreaPhi()));
630  pj[2] = jet->Pz() - jet->Area()*(rho*TMath::SinH(jet->AreaEta()));
631 
632  return Z(p,pj);
633 }
634 //_______________________________________________________________________________
635 
636 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
637 
638  Double_t p[3],pj[3];
639  Bool_t okpp=part->PxPyPz(p);
640  Bool_t okpj=jet->PxPyPz(pj);
641 
642  if(!okpp || !okpj)
643  {
644  printf("Problems getting momenta\n");
645  return -999;
646  }
647 
648  return Z(p,pj);
649 }
650 //_______________________________________________________________________________
651 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
652 
653  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
654  Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
655  return z;
656 }
657 
658 
659 //_______________________________________________________________________________
660 Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
661 
662  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
663  Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
664  return z;
665 }
666 
667 //_______________________________________________________________________________
668 
670 
671  // Statistics
672  Int_t nbins=8;
673  if(fUseMCInfo) nbins+=2;
674 
675  fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
676  fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
677  fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
678  fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
679  fhstat->GetXaxis()->SetBinLabel(4,"N jets");
680  fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
681  fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
682  fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
683  fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
684  if(fUseMCInfo) {
685  fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
686  fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
687  fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
688  fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
689  }
690  fhstat->SetNdivisions(1);
691  fOutput->Add(fhstat);
692 
693  fhCentDjet=new TH1F("hCentDjet","Centrality",100,0,100);
694  fOutput->Add(fhCentDjet);
695 
696  const Int_t nbinsmass=300;
697  const Int_t nbinspttrack=500;
698  Int_t nbinsptjet=200;
699  const Int_t nbinsptD=100;
700  const Int_t nbinsphi=200;
701  const Int_t nbinseta=100;
702 
703  //binning for THnSparse
704  const Int_t nbinsSpsmass=60;
705  const Int_t nbinsSpsptjet=200;
706  const Int_t nbinsSpsptD=50;
707  const Int_t nbinsSpsz=144;
708 
709  const Float_t pttracklims[2]={0.,200.};
710  Float_t ptjetlims[2]={-50.,150.};
711  const Float_t ptDlims[2]={0.,50.};
712  const Float_t zlims[2]={0.,1.2};
713  const Float_t philims[2]={0.,6.3};
714  const Float_t etalims[2]={-1.5,1.5};
715 
716  // jet related fistograms
717 
718  fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
719  fhPhiJetTrks->Sumw2();
720  fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
721  fhEtaJetTrks->Sumw2();
722  fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
723  fhPtJetTrks->Sumw2();
724 
725  fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
726  fhPhiJet->Sumw2();
727  fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
728  fhEtaJet->Sumw2();
729  fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
730  fhPtJet->Sumw2();
731 
732  fOutput->Add(fhPhiJetTrks);
733  fOutput->Add(fhEtaJetTrks);
734  fOutput->Add(fhPtJetTrks);
735  fOutput->Add(fhPhiJet);
736  fOutput->Add(fhEtaJet);
737  fOutput->Add(fhPtJet);
738 
740  {
741  if(fAnalyseDBkg){
742  fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
743  fhDiffSideBand->SetStats(kTRUE);
744  fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
745  fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
746  fhDiffSideBand->Sumw2();
747  fOutput->Add(fhDiffSideBand);
748  }
749 
750  fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
751  fhPtPion->SetStats(kTRUE);
752  fhPtPion->GetXaxis()->SetTitle("GeV/c");
753  fhPtPion->GetYaxis()->SetTitle("Entries");
754  fhPtPion->Sumw2();
755  fOutput->Add(fhPtPion);
756 
757  }
758 
759 
760  fhInvMassptD = new TH2F("hInvMassptD","D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
761  fhInvMassptD->SetStats(kTRUE);
762  fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
763  fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
764  fhInvMassptD->Sumw2();
765 
766  fOutput->Add(fhInvMassptD);
767 
768  if(fUseMCInfo){
769  fhInvMassptDbg = new TH2F("hInvMassptDbg","Background D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
770  fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
771  fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
772  fhInvMassptDbg->Sumw2();
773  fOutput->Add(fhInvMassptDbg);
774 
775  }
776 
777  fhsDphiz=0x0; //definition below according to the switches
778 
779  if(fUseMCInfo){
780  AliInfo("Creating a 7 axes container (MB background candidates)");
781  const Int_t nAxis=7;
782  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
783  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
784  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
785  fNAxesBigSparse=nAxis;
786  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
787 
788  } else{
789  AliInfo("Creating a 5 axes container");
790  const Int_t nAxis=5;
791  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,63};
792  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,0};
793  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,6.3};
794  fNAxesBigSparse=nAxis;
795 
796  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., Jet phi", nAxis, nbinsSparse, minSparse, maxSparse);
797  }
798 
799  if(!fhsDphiz) AliFatal("No THnSparse created");
800  fhsDphiz->Sumw2();
801 
802  fOutput->Add(fhsDphiz);
803 
805  {
806  const Int_t RMnbins[6] = {nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsz,nbinsSpsptjet,nbinsSpsptD};
807  const Double_t RMmin[6]={zlims[0],ptjetlims[0],ptDlims[0],zlims[0],ptjetlims[0],ptDlims[0]};
808  const Double_t RMmax[6]={zlims[1],ptjetlims[1],ptDlims[1],zlims[1],ptjetlims[1],ptDlims[1]};
809  fResponseMatrix = new THnSparseF("ResponseMatrix","z, jet pt, D pt: Rec and Gen",6,RMnbins,RMmin,RMmax);
810  fResponseMatrix->Sumw2();
811  fOutput->Add(fResponseMatrix);
812  }
813 
814  PostData(1, fOutput);
815 
816  return kTRUE;
817 }
818 
819 //_______________________________________________________________________________
820 
821 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent){
822 
823 
824  Double_t masses[2]={0.,0.};
825  Int_t pdgdaughtersD0[2]={211,321};//pi,K
826  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
827 
828  masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
829  masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
830 
831  Double_t *point=0x0;
832 
833  if(!fAnalyseDBkg)
834  {
835  point=new Double_t[5];
836  point[0]=z;
837  point[1]=ptj;
838  point[2]=ptD;
839  point[3]=masses[0];
840  point[4]=phij;
841 
842  }
843  else
844  {
845  point=new Double_t[7];
846  point[0]=z;
847  point[1]=ptj;
848  point[2]=ptD;
849  point[3]=masses[0];
850  point[4]=static_cast<Double_t>(IsBkg ? 1 : 0);
851  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
852  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
853 
854  }
855  Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
856  if(isselected==1 || isselected==3)
857  {
858 
859  fhInvMassptD->Fill(masses[0],ptD);
860  fhsDphiz->Fill(point,1.);
861 
862  }
863  if(isselected>=2)
864  {
865  fhInvMassptD->Fill(masses[1],ptD);
866  point[3]=masses[1];
867  fhsDphiz->Fill(point,1.);
868 
869  }
870  delete[] point;
871 }
872 
873 //_______________________________________________________________________________
874 
875 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
876 
877  AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
878  Double_t deltamass= dstar->DeltaInvMass();
879 
880 
881  fhPtPion->Fill(softpi->Pt());
882 
883  Double_t *point=0x0;
884  if(!fAnalyseDBkg)
885  {
886  point=new Double_t[5];
887  point[0]=z;
888  point[1]=ptj;
889  point[2]=ptD;
890  point[3]=deltamass;
891  point[4]=phij;
892  }
893  else
894  {
895  point=new Double_t[7];
896  point[0]=z;
897  point[1]=ptj;
898  point[2]=ptD;
899  point[3]=deltamass;
900  point[4]=static_cast<Double_t>(IsBkg ? 1 : 0);
901  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
902  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
903  }
904 
905  if(!point){
906  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
907  return;
908  }
909 
910 
911  fhInvMassptD->Fill(deltamass,ptD);
912  fhsDphiz->Fill(point,1.);
913  delete[] point;
914 }
915 
916 //_______________________________________________________________________________
917 
918 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t z,Double_t ptD,Double_t ptjet, Double_t phij, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
919 
920  Double_t pdgmass=0;
921  if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
922  if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
923 
924  Double_t *point=0x0;
925 
926  if(fNAxesBigSparse==5){
927  point=new Double_t[5];
928  point[0]=z;
929  point[1]=ptjet;
930  point[2]=ptD;
931  point[3]=pdgmass;
932  point[4]=phij;
933  }
934  if(fNAxesBigSparse==7){
935  point=new Double_t[7];
936  point[0]=z;
937  point[1]=ptjet;
938  point[2]=ptD;
939  point[3]=pdgmass;
940  point[4]=1;
941  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
942  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
943  }
944 
945  if(!point){
946  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
947  return;
948  }
949 
950 
951  point[3]=pdgmass;
952  fhsDphiz->Fill(point,1.);
953 
954  delete[] point;
955 }
956 
957 //_______________________________________________________________________________
958 
959 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2, Double_t rho) const {
960  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
961  //It recalculates the eta-phi values if it was asked for background subtraction of the jets
962  if(!p1 || !p2) return -1;
963 
964  Double_t phi1=p1->Phi(),eta1=p1->Eta();
965 
966  if (rho>0)
967  {
968  Double_t pj[3];
969  Bool_t okpj=p1->PxPyPz(pj);
970  if(!okpj){
971  printf("Problems getting momenta\n");
972  return -999;
973  }
974  pj[0] = p1->Px() - p1->Area()*(rho*TMath::Cos(p1->AreaPhi()));
975  pj[1] = p1->Py() - p1->Area()*(rho*TMath::Sin(p1->AreaPhi()));
976  pj[2] = p1->Pz() - p1->Area()*(rho*TMath::SinH(p1->AreaEta()));
977  //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
978  //[0,pi] if py > 0 and
979  //[pi,2pi] if py < 0
980  phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
981  if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
982  eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
983  }
984 
985  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
986 
987  Double_t dPhi=TMath::Abs(phi1-phi2);
988  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
989 
990 
991  Double_t dEta=eta1-eta2;
992  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
993  return deltaR;
994 
995 }
996 //_______________________________________________________________________________
998  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
999  if(!p1 || !p2) return -1;
1000 
1001  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1002 
1003  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1004 
1005  Double_t dPhi=TMath::Abs(phi1-phi2);
1006  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1007 
1008 
1009  Double_t dEta=eta1-eta2;
1010  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1011  return deltaR;
1012 
1013 }
1014 
1015 //_______________________________________________________________________________
1016 
1018 
1019  Double_t ptD=candDstar->Pt();
1020  Int_t bin = fCuts->PtBin(ptD);
1021  if (bin < 0){
1022  // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1023  bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1024  return -1; // out of bounds
1025  }
1026 
1027  Double_t invM=candDstar->InvMassD0();
1028  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1029 
1030  Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1031  Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1032 
1033  if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1034  else return 0;
1035 
1036 }
1037 
1038 //_______________________________________________________________________________
1039 
1041  //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1042 
1043  Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1044  Bool_t binEMCal=kTRUE;
1045  Double_t phi=vpart->Phi(), eta=vpart->Eta();
1046  if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1047  if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1048  return binEMCal;
1049 
1050 
1051 }
Int_t pdg
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:217
Double_t Area() const
Definition: AliEmcalJet.h:97
Double_t GetRhoVal() const
const TString & GetRhoName() const
Double_t DeltaInvMass() const
Jet is tagged to contain a D* meson.
Definition: AliEmcalJet.h:52
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:88
Double_t Z(AliVParticle *part, AliEmcalJet *jet, Double_t rho) const
Double_t Py() const
Definition: AliEmcalJet.h:74
Double_t Phi() const
Definition: AliEmcalJet.h:84
Double_t mass
Int_t GetNParticles() const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
AliVTrack * GetTrack() const
TList * fOutput
!output list
void AngularCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
AliVParticle * GetFlavourTrack(Int_t i=0) const
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:106
Double_t Px() const
Definition: AliEmcalJet.h:73
AliParticleContainer * GetParticleContainer() const
Double_t AreaEta() const
Definition: AliEmcalJet.h:99
void FillDJetHistograms(AliEmcalJet *jet, Double_t rho, Bool_t IsBkg, AliAODEvent *aodEvent)
void AddFlavourTrack(AliVParticle *hftrack)
virtual AliVParticle * GetParticle(Int_t i=-1) const
AliAODTrack * GetBachelor() const
Double_t fCent
!event centrality
TClonesArray * GetArray() const
Float_t CheckDeltaR(AliEmcalJet *p1, AliVParticle *p2) const
ClassImp(AliAnalysisTaskFlavourJetCorrelations) AliAnalysisTaskFlavourJetCorrelations
Double_t AreaPhi() const
Definition: AliEmcalJet.h:100
Double_t Pt() const
Definition: AliEmcalJet.h:76
AliEmcalJet * GetNextJet()
Bool_t IsEventSelected(AliVEvent *event)
Jet is tagged to contain a D0 meson.
Definition: AliEmcalJet.h:53
Float_t GetJetRadius() const
void FillHistogramsMCGenDJetCorr(Double_t z, Double_t ptD, Double_t ptjet, Double_t phij, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
void ConstituentCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:127
void FillHistogramsD0JetCorr(AliAODRecoDecayHF *candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent *aodEvent)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:274
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:78
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
Double_t Pz() const
Definition: AliEmcalJet.h:75
void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF *dstar, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
Double_t InvMassD0() const
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:233
const Int_t nbins
Int_t PtBin(Double_t pt) const
void ResetCurrentID(Int_t i=-1)
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...