1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Analysis Taks for reconstructed particle correlation
18 // (first implementation done for D mesons) with jets
19 // (use the so called Emcal framework)
21 //-----------------------------------------------------------------------
23 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
24 // A. Grelli (Utrecht University) a.grelli@uu.nl
25 // X. Zhang (LBNL) XMZhang@lbl.gov
26 //-----------------------------------------------------------------------
28 #include <TDatabasePDG.h>
29 #include <TParticle.h>
32 #include <THnSparse.h>
34 #include "AliAnalysisTaskFlavourJetCorrelations.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODHandler.h"
37 #include "AliAnalysisManager.h"
39 #include "AliEmcalJet.h"
40 #include "AliJetContainer.h"
41 #include "AliAODRecoDecay.h"
42 #include "AliAODRecoCascadeHF.h"
43 #include "AliAODRecoDecayHF2Prong.h"
44 #include "AliESDtrack.h"
45 #include "AliAODMCParticle.h"
46 #include "AliPicoTrack.h"
47 #include "AliRDHFCutsD0toKpi.h"
48 #include "AliRDHFCutsDStartoKpipi.h"
50 ClassImp(AliAnalysisTaskFlavourJetCorrelations)
53 //_______________________________________________________________________________
55 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
56 AliAnalysisTaskEmcalJet("",kFALSE),
70 fLeadingJetOnly(kFALSE),
78 //SetMakeGeneralHistograms(kTRUE);
82 //_______________________________________________________________________________
84 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
85 AliAnalysisTaskEmcalJet(name,kFALSE),
99 fLeadingJetOnly(kFALSE),
105 // Constructor. Initialization of Inputs and Outputs
108 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
110 fCandidateType=candtype;
111 const Int_t nptbins=fCuts->GetNPtBins();
112 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};
114 switch(fCandidateType){
118 fPDGdaughters[0]=211;//pi
119 fPDGdaughters[1]=321;//K
120 fPDGdaughters[2]=0; //empty
121 fPDGdaughters[3]=0; //empty
122 fBranchName="D0toKpi";
128 fPDGdaughters[1]=211;//pi soft
129 fPDGdaughters[0]=421;//D0
130 fPDGdaughters[2]=211;//pi fromD0
131 fPDGdaughters[3]=321; // K from D0
133 fCandArrName="Dstar";
136 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
138 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
142 printf("%d not accepted!!\n",fCandidateType);
146 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
147 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
149 DefineInput(1, TClonesArray::Class());
150 DefineInput(2, TClonesArray::Class());
152 DefineOutput(1,TList::Class()); // histos
153 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
157 //_______________________________________________________________________________
159 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
164 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
171 //_______________________________________________________________________________
173 void AliAnalysisTaskFlavourJetCorrelations::Init(){
178 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
179 switch(fCandidateType){
182 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
183 copyfCuts->SetName("AnalysisCutsDzero");
185 PostData(2,copyfCuts);
190 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
191 copyfCuts->SetName("AnalysisCutsDStar");
193 PostData(2,copyfCuts);
203 //_______________________________________________________________________________
205 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
207 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
208 fmyOutput = new TList();
209 fmyOutput->SetOwner();
210 fmyOutput->SetName("pippo");
212 DefineHistoForAnalysis();
213 PostData(1,fmyOutput);
218 //_______________________________________________________________________________
220 void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)
224 AliAnalysisTaskEmcalJet::ExecOnce();
228 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
230 TClonesArray *arrayDStartoD0pi=0;
231 TClonesArray *trackArr = 0;
233 if (!aodEvent && AODEvent() && IsStandardAOD()) {
235 // In case there is an AOD handler writing a standard AOD, use the AOD
236 // event in memory rather than the input (ESD) event.
237 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
239 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
240 // have to taken from the AOD event hold by the AliAODExtension
241 AliAODHandler* aodHandler = (AliAODHandler*)
242 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
243 if(aodHandler->GetExtensions()) {
244 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
245 AliAODEvent *aodFromExt = ext->GetAOD();
246 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
249 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
252 if (!arrayDStartoD0pi) {
253 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
255 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
257 TClonesArray* mcArray = 0x0;
259 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
261 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
267 trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
268 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
269 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
270 fJetRadius=GetJetContainer(0)->GetJetRadius();
273 AliInfo(Form("Could not find the track array\n"));
278 if (!fCandidateArray) fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
279 if (fCandidateType==1 && !fSideBandArray) fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
281 Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
284 TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
285 TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
286 TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
287 TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
288 TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
289 TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
290 TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
291 TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
292 TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
293 TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
294 TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
295 TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
296 TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
297 TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
298 TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
299 TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
300 THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
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;
309 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
310 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
311 if(!iseventselected) return;
315 //retrieve charm candidates selected
316 Int_t candidates = fCandidateArray->GetEntriesFast();
320 Int_t njets=GetJetContainer()->GetNJets();
322 hstat->Fill(6, candidates);
323 hNDPerEvNoJet->Fill(candidates);
324 for(Int_t iD=0;iD<candidates;iD++){
325 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
327 hptDPerEvNoJet->Fill(cand->Pt());
334 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
336 for(Int_t ic = 0; ic < candidates; ic++) {
339 AliAODRecoDecayHF* charm=0x0;
340 AliAODRecoCascadeHF* dstar=0x0;
343 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
346 if(fCandidateType==kDstartoKpipi){
347 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
353 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
355 if(fCandidateType==kDstartoKpipi) {
356 Double_t deltamass= dstar->DeltaInvMass();
357 candsparse[3]=deltamass;
358 hnspDstandalone->Fill(candsparse);
360 if(fCandidateType==kD0toKpi){
361 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
362 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
365 Int_t pdgdaughtersD0[2]={211,321};//pi,K
366 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
368 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
369 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
370 if(isselected==1 || isselected==3) {
371 candsparse[3]=masses[0];
372 hnspDstandalone->Fill(candsparse);
375 candsparse[3]=masses[1];
376 hnspDstandalone->Fill(candsparse);
382 // we start with jets
387 Double_t leadingJet =0;
389 Int_t ntrarr=trackArr->GetEntriesFast();
390 hNtrArr->Fill(ntrarr);
392 for(Int_t i=0;i<ntrarr;i++){
393 AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
395 hPtJetTrks->Fill(jtrack->Pt());
396 hPhiJetTrks->Fill(jtrack->Phi());
397 hEtaJetTrks->Fill(jtrack->Eta());
398 hEjetTrks->Fill(jtrack->E());
402 //option to use only the leading jet
404 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
405 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
407 if(ptjet>leadingJet ) leadingJet = ptjet;
413 //loop over jets and consider only the leading jet in the event
414 for (Int_t iJets = 0; iJets<njets; iJets++) {
415 //Printf("Jet N %d",iJets);
416 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
417 //Printf("Corr task Accept Jet");
418 if(!AcceptJet(jet)) {
429 // choose the leading jet
430 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
431 //used for normalization
434 // fill energy, eta and phi of the jet
436 hPhiJet ->Fill(phiJet);
437 hEtaJet ->Fill(etaJet);
438 hPtJet ->Fill(ptjet);
440 //loop on jet particles
441 Int_t ntrjet= jet->GetNumberOfTracks();
442 for(Int_t itrk=0;itrk<ntrjet;itrk++){
444 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
445 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
447 }//end loop on jet tracks
451 hPtJetPerEvNoD->Fill(jet->Pt());
454 //Printf("N candidates %d ", candidates);
455 for(Int_t ic = 0; ic < candidates; ic++) {
458 AliVParticle* charm=0x0;
459 charm=(AliVParticle*)fCandidateArray->At(ic);
463 FlagFlavour(charm, jet);
464 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
466 FillHistogramsRecoJetCorr(charm, jet);
469 //retrieve side band background candidates for Dstar
470 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
472 for(Int_t ib=0;ib<nsbcand;ib++){
473 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
474 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
475 if(!sbcand) continue;
476 SideBandBackground(sbcand,jet);
479 AliAODRecoDecayHF* charmbg = 0x0;
480 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
481 if(!charmbg) continue;
482 MCBackground(charmbg,jet);
487 hNJetPerEv->Fill(cntjet);
488 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
489 PostData(1,fmyOutput);
493 //_______________________________________________________________________________
495 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
497 // The Terminate() function is the last function to be called during
498 // a query. It always runs on the client, it can be used to present
499 // the results graphically or save the results to file.
501 Info("Terminate"," terminate");
502 AliAnalysisTaskSE::Terminate();
504 fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
506 printf("ERROR: fmyOutput not available\n");
511 //_______________________________________________________________________________
513 void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
515 Int_t abspdg=TMath::Abs(pdg);
517 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
518 // compute the Delta mass for the D*
519 if(fCandidateType==kDstartoKpipi){
521 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
525 fMinMass = mass-range;
526 fMaxMass = mass+range;
528 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
529 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
532 //_______________________________________________________________________________
534 void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
541 printf("Error! Lower limit larger than upper limit!\n");
542 fMinMass = uplimit - uplimit*0.2;
547 //_______________________________________________________________________________
549 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
551 AliInfo("Maximum number of bins allowed is 30!");
554 if(!width) return kFALSE;
555 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
559 //_______________________________________________________________________________
561 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
563 printf("Particle or jet do not exist!\n");
567 Bool_t okpp=part->PxPyPz(p);
568 Bool_t okpj=jet->PxPyPz(pj);
570 printf("Problems getting momenta\n");
573 Double_t pjet=jet->P();
574 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
578 //_______________________________________________________________________________
580 Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
583 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
584 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
585 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
586 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
587 hstat->GetXaxis()->SetBinLabel(4,"N jets");
588 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
589 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
590 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
591 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
592 hstat->SetNdivisions(1);
593 fmyOutput->Add(hstat);
595 const Int_t nbinsmass=200;
596 const Int_t nbinsptjet=500;
597 const Int_t nbinsptD=100;
598 const Int_t nbinsz=100;
599 const Int_t nbinsphi=200;
600 const Int_t nbinseta=100;
602 const Float_t ptjetlims[2]={0.,200.};
603 const Float_t ptDlims[2]={0.,50.};
604 const Float_t zlims[2]={0.,1.2};
605 const Float_t philims[2]={0.,6.3};
606 const Float_t etalims[2]={-1.5,1.5};
608 if(fCandidateType==kDstartoKpipi)
611 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
612 hDiffSideBand->SetStats(kTRUE);
613 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
614 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
615 hDiffSideBand->Sumw2();
616 fmyOutput->Add(hDiffSideBand);
619 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
620 hPtPion->SetStats(kTRUE);
621 hPtPion->GetXaxis()->SetTitle("GeV/c");
622 hPtPion->GetYaxis()->SetTitle("Entries");
624 fmyOutput->Add(hPtPion);
627 // D related histograms
628 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
629 hNDPerEvNoJet->Sumw2();
630 fmyOutput->Add(hNDPerEvNoJet);
632 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
633 hptDPerEvNoJet->Sumw2();
634 fmyOutput->Add(hptDPerEvNoJet);
636 const Int_t nAxisD=4;
637 const Int_t nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
638 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
639 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
640 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
641 hsDstandalone->Sumw2();
643 fmyOutput->Add(hsDstandalone);
645 // jet related fistograms
647 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
649 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
650 hPhiJetTrks->Sumw2();
651 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
652 hEtaJetTrks->Sumw2();
653 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
656 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
658 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
660 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
662 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
665 TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
666 hPtJetWithD->Sumw2();
667 //for the MC this histogram is filled with the real background
668 TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
669 hPtJetWithDsb->Sumw2();
671 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
672 hdeltaRJetTracks->Sumw2();
674 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
677 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
680 TH1F *hNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
681 hNJetPerEvNoD->Sumw2();
683 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
684 hPtJetPerEvNoD->Sumw2();
686 fmyOutput->Add(hEjetTrks);
687 fmyOutput->Add(hPhiJetTrks);
688 fmyOutput->Add(hEtaJetTrks);
689 fmyOutput->Add(hPtJetTrks);
690 fmyOutput->Add(hEjet);
691 fmyOutput->Add(hPhiJet);
692 fmyOutput->Add(hEtaJet);
693 fmyOutput->Add(hPtJet);
694 fmyOutput->Add(hPtJetWithD);
695 fmyOutput->Add(hPtJetWithDsb);
696 fmyOutput->Add(hdeltaRJetTracks);
697 fmyOutput->Add(hNtrArr);
698 fmyOutput->Add(hNJetPerEv);
699 fmyOutput->Add(hNJetPerEvNoD);
700 fmyOutput->Add(hPtJetPerEvNoD);
702 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
704 fmyOutput->Add(hDeltaRD);
706 //background (side bands for the Dstar and like sign for D0)
707 fJetRadius=GetJetContainer(0)->GetJetRadius();
708 TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
709 hInvMassptD->SetStats(kTRUE);
710 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
711 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
712 hInvMassptD->Sumw2();
714 fmyOutput->Add(hInvMassptD);
717 TH2F* hInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
718 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
719 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
720 hInvMassptDbg->Sumw2();
721 fmyOutput->Add(hInvMassptDbg);
724 Double_t pi=TMath::Pi(), philims2[2];
725 philims2[0]=-pi*1./2.;
726 philims2[1]=pi*3./2.;
728 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2};
729 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5};
730 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5};
731 THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass and IsSB", nAxis, nbinsSparse, minSparse, maxSparse);
734 fmyOutput->Add(hsDphiz);
736 PostData(1, fmyOutput);
741 //_______________________________________________________________________________
743 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
745 Double_t ptD=candidate->Pt();
746 Double_t ptjet=jet->Pt();
747 Double_t deltaR=DeltaR(candidate,jet);
748 Double_t deltaphi = jet->Phi()-candidate->Phi();
749 if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
750 if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
751 Double_t z=Z(candidate,jet);
753 TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
754 hDeltaRD->Fill(deltaR);
756 if(fCandidateType==kD0toKpi) {
757 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
758 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
762 if(fCandidateType==kDstartoKpipi) {
763 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
764 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
767 } else{ //generated level
768 //AliInfo("Non reco");
769 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
773 //_______________________________________________________________________________
775 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
778 Double_t masses[2]={0.,0.};
779 Int_t pdgdaughtersD0[2]={211,321};//pi,K
780 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
782 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
783 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
785 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
786 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
787 Double_t point[6]={z,dPhi,ptj,ptD,masses[0],0};
789 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
790 if(isselected==1 || isselected==3) {
792 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
794 FillMassHistograms(masses[0], ptD, deltaR);
795 hsDphiz->Fill(point,1.);
798 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
800 FillMassHistograms(masses[1], ptD, deltaR);
802 hsDphiz->Fill(point,1.);
807 //_______________________________________________________________________________
809 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
810 //dPhi and z not used at the moment,but will be (re)added
812 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
813 Double_t deltamass= dstar->DeltaInvMass();
814 //Double_t massD0= dstar->InvMassD0();
817 TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
818 hPtPion->Fill(softpi->Pt());
820 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
821 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
822 Int_t isSB=IsDzeroSideBand(dstar);
824 Double_t point[6]={z,dPhi,ptj,ptD,deltamass,isSB};
826 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
828 FillMassHistograms(deltamass, ptD, deltaR);
829 hsDphiz->Fill(point,1.);
833 //_______________________________________________________________________________
835 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
838 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
839 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
840 Double_t point[6]={z,dPhi,ptjet,ptD,pdgmass,0};
842 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
843 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
846 if(deltaR<fJetRadius) {
847 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
848 hsDphiz->Fill(point,1.);
853 //_______________________________________________________________________________
855 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
857 if(deltaR<fJetRadius) {
858 TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
859 hInvMassptD->Fill(mass,ptD);
863 //________________________________________________________________________________
865 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
866 Double_t deltaR=DeltaR(charm, jet);
867 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
868 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
869 if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
875 //_______________________________________________________________________________
877 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
879 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
880 // (expected detector resolution) on the left and right frm the D0 mass. Each band
881 // has a width of ~5 sigmas. Two band needed for opening angle considerations
882 TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
883 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
885 Double_t deltaM=candDstar->DeltaInvMass();
886 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
887 Double_t ptD=candDstar->Pt();
888 Double_t ptjet=jet->Pt();
889 Double_t dPhi=jet->Phi()-candDstar->Phi();
890 Double_t deltaR=DeltaR(candDstar,jet);
891 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
892 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
894 Int_t isSideBand=IsDzeroSideBand(candDstar);
895 //fill the background histogram with the side bands or when looking at MC with the real background
897 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
898 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
900 if(deltaR<fJetRadius){ // evaluate in the near side
901 //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
902 hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
907 //_______________________________________________________________________________
909 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
911 //need mass, deltaR, pt jet, ptD
912 //two cases: D0 and Dstar
914 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
915 if(!isselected) return;
917 Double_t ptD=candbg->Pt();
918 Double_t ptjet=jet->Pt();
919 Double_t deltaR=DeltaR(candbg,jet);
921 TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
922 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
925 if(fCandidateType==kDstartoKpipi){
926 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
927 Double_t deltaM=dstarbg->DeltaInvMass();
928 hInvMassptDbg->Fill(deltaM,ptD);
929 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
932 if(fCandidateType==kD0toKpi){
933 Double_t masses[2]={0.,0.};
934 Int_t pdgdaughtersD0[2]={211,321};//pi,K
935 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
937 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
938 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
941 if(isselected==1 || isselected==3) {
942 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
943 hInvMassptDbg->Fill(masses[0],ptD);
946 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
947 hInvMassptDbg->Fill(masses[1],ptD);
954 //_______________________________________________________________________________
956 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
957 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
959 if(!p1 || !p2) return -1;
960 Double_t phi1=p1->Phi(),eta1=p1->Eta();
961 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
963 Double_t dPhi=phi1-phi2;
964 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
965 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
967 Double_t dEta=eta1-eta2;
968 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
973 //_______________________________________________________________________________
975 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
977 Double_t ptD=candDstar->Pt();
978 Int_t bin = fCuts->PtBin(ptD);
980 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
981 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
982 return -1; // out of bounds
985 Double_t invM=candDstar->InvMassD0();
986 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
988 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
989 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
991 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;