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