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("",kTRUE),
69 fLeadingJetOnly(kFALSE),
83 //SetMakeGeneralHistograms(kTRUE);
87 //_______________________________________________________________________________
89 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, Bool_t jetOnly) :
90 AliAnalysisTaskEmcalJet(name,kTRUE),
103 fLeadingJetOnly(kFALSE),
107 fJetOnlyMode(jetOnly),
115 // Constructor. Initialization of Inputs and Outputs
118 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
120 fCandidateType=candtype;
121 const Int_t nptbins=fCuts->GetNPtBins();
122 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};
124 switch(fCandidateType){
128 fPDGdaughters[0]=211;//pi
129 fPDGdaughters[1]=321;//K
130 fPDGdaughters[2]=0; //empty
131 fPDGdaughters[3]=0; //empty
132 fBranchName="D0toKpi";
138 fPDGdaughters[1]=211;//pi soft
139 fPDGdaughters[0]=421;//D0
140 fPDGdaughters[2]=211;//pi fromD0
141 fPDGdaughters[3]=321; // K from D0
143 fCandArrName="Dstar";
146 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
148 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
152 printf("%d not accepted!!\n",fCandidateType);
156 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
157 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
160 DefineOutput(1,TList::Class()); //histos with jet info
161 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
164 DefineInput(1, TClonesArray::Class());
165 DefineInput(2, TClonesArray::Class());
167 DefineOutput(1,TList::Class()); // histos
168 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
172 //_______________________________________________________________________________
174 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
179 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
185 //_______________________________________________________________________________
187 void AliAnalysisTaskFlavourJetCorrelations::Init(){
192 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
194 switch(fCandidateType){
197 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
198 copyfCuts->SetName("AnalysisCutsDzero");
200 PostData(2,copyfCuts);
205 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
206 copyfCuts->SetName("AnalysisCutsDStar");
208 PostData(2,copyfCuts);
218 //_______________________________________________________________________________
220 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
222 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
223 AliAnalysisTaskEmcal::UserCreateOutputObjects();
225 // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
226 DefineHistoForAnalysis();
232 //_______________________________________________________________________________
234 Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
236 // user exec from AliAnalysisTaskEmcal is used
239 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
241 TClonesArray *arrayDStartoD0pi=0;
243 if (!aodEvent && AODEvent() && IsStandardAOD()) {
245 // In case there is an AOD handler writing a standard AOD, use the AOD
246 // event in memory rather than the input (ESD) event.
247 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
249 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
250 // have to taken from the AOD event hold by the AliAODExtension
251 AliAODHandler* aodHandler = (AliAODHandler*)
252 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
253 if(aodHandler->GetExtensions()) {
254 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
255 AliAODEvent *aodFromExt = ext->GetAOD();
256 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
259 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
262 if (!arrayDStartoD0pi) {
263 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
265 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
267 TClonesArray* mcArray = 0x0;
269 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
271 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
277 fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
278 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
279 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
280 fJetRadius=GetJetContainer(0)->GetJetRadius();
283 AliInfo(Form("Could not find the track array\n"));
288 fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
289 if (!fCandidateArray) return kFALSE;
290 if (fCandidateType==1) {
291 fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
292 if (!fSideBandArray) return kFALSE;
294 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
297 TH1I* hstat=(TH1I*)fOutput->FindObject("hstat");
298 TH1F* hPtJetTrks=(TH1F*)fOutput->FindObject("hPtJetTrks");
299 TH1F* hPhiJetTrks=(TH1F*)fOutput->FindObject("hPhiJetTrks");
300 TH1F* hEtaJetTrks=(TH1F*)fOutput->FindObject("hEtaJetTrks");
301 TH1F* hEjetTrks=(TH1F*)fOutput->FindObject("hEjetTrks");
302 TH1F* hPtJet=(TH1F*)fOutput->FindObject("hPtJet");
303 TH1F* hPhiJet=(TH1F*)fOutput->FindObject("hPhiJet");
304 TH1F* hEtaJet=(TH1F*)fOutput->FindObject("hEtaJet");
305 TH1F* hEjet=(TH1F*)fOutput->FindObject("hEjet");
306 TH1F* hNtrArr=(TH1F*)fOutput->FindObject("hNtrArr");
307 TH1F* hNJetPerEv=(TH1F*)fOutput->FindObject("hNJetPerEv");
308 TH1F* hdeltaRJetTracks=(TH1F*)fOutput->FindObject("hdeltaRJetTracks");
309 TH1F* hNDPerEvNoJet=(TH1F*)fOutput->FindObject("hNDPerEvNoJet");
310 TH1F* hptDPerEvNoJet=(TH1F*)fOutput->FindObject("hptDPerEvNoJet");
311 TH1F* hNJetPerEvNoD=(TH1F*)fOutput->FindObject("hNJetPerEvNoD");
312 TH1F* hPtJetPerEvNoD=(TH1F*)fOutput->FindObject("hPtJetPerEvNoD");
313 THnSparseF* hnspDstandalone=(THnSparseF*)fOutput->FindObject("hsDstandalone");
314 THnSparseF* hsJet=(THnSparseF*)fOutput->FindObject("hsJet");
318 // fix for temporary bug in ESDfilter
319 // the AODs with null vertex pointer didn't pass the PhysSel
320 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
323 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
324 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
325 if(!iseventselected) return kFALSE;
329 //retrieve charm candidates selected
330 Int_t candidates = 0;
331 Int_t njets=GetJetContainer()->GetNJets();
334 candidates = fCandidateArray->GetEntriesFast();
338 hstat->Fill(6, candidates);
339 hNDPerEvNoJet->Fill(candidates);
340 for(Int_t iD=0;iD<candidates;iD++){
341 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
343 hptDPerEvNoJet->Fill(cand->Pt());
350 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
352 for(Int_t ic = 0; ic < candidates; ic++) {
355 AliAODRecoDecayHF* charm=0x0;
356 AliAODRecoCascadeHF* dstar=0x0;
359 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
362 if(fCandidateType==kDstartoKpipi){
363 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
369 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
371 if(fCandidateType==kDstartoKpipi) {
373 Double_t deltamass= dstar->DeltaInvMass();
374 candsparse[3]=deltamass;
375 } else candsparse[3] = 0.145; //for generated
376 hnspDstandalone->Fill(candsparse);
378 if(fCandidateType==kD0toKpi){
380 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
381 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
382 //not a further selection,just to get the value of isselected for the D0
384 Int_t pdgdaughtersD0[2]={211,321};//pi,K
385 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
387 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
388 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
389 if(isselected==1 || isselected==3) {
390 candsparse[3]=masses[0];
391 hnspDstandalone->Fill(candsparse);
394 candsparse[3]=masses[1];
395 hnspDstandalone->Fill(candsparse);
399 Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
400 candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
401 hnspDstandalone->Fill(candsparse);
407 // we start with jets
412 Double_t leadingJet =0;
415 Int_t ntrarr=fTrackArr->GetEntriesFast();
416 hNtrArr->Fill(ntrarr);
418 for(Int_t i=0;i<ntrarr;i++){
419 AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
421 hPtJetTrks->Fill(jtrack->Pt());
422 hPhiJetTrks->Fill(jtrack->Phi());
423 hEtaJetTrks->Fill(jtrack->Eta());
424 hEjetTrks->Fill(jtrack->E());
428 //option to use only the leading jet
430 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
431 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
433 if(ptjet>leadingJet ) leadingJet = ptjet;
439 //loop over jets and consider only the leading jet in the event
440 for (Int_t iJets = 0; iJets<njets; iJets++) {
441 //Printf("Jet N %d",iJets);
442 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
443 //Printf("Corr task Accept Jet");
444 if(!AcceptJet(jet)) {
454 Double_t origPtJet=fPtJet;
459 pointJ[4] = jet->GetNumberOfConstituents();
460 pointJ[5] = jet->Area();
462 // choose the leading jet
463 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
464 //used for normalization
467 // fill energy, eta and phi of the jet
469 hPhiJet ->Fill(phiJet);
470 hEtaJet ->Fill(etaJet);
471 hPtJet ->Fill(fPtJet);
472 if(fJetOnlyMode) hsJet->Fill(pointJ,1);
473 //loop on jet particles
474 Int_t ntrjet= jet->GetNumberOfTracks();
475 for(Int_t itrk=0;itrk<ntrjet;itrk++){
477 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
478 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
480 }//end loop on jet tracks
484 hPtJetPerEvNoD->Fill(fPtJet);
487 //Printf("N candidates %d ", candidates);
488 for(Int_t ic = 0; ic < candidates; ic++) {
491 AliVParticle* charm=0x0;
492 charm=(AliVParticle*)fCandidateArray->At(ic);
494 AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
495 fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
496 if (fIsDInJet) FlagFlavour(jet);
497 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
499 //Note: the z component of the jet momentum comes from the eta-phi direction of the jet particles, it is not calculated from the z component of the tracks since, as default, the scheme used for jet reco is the pt-scheme which sums the scalar component, not the vectors. Addind the D daughter momentum component by componet as done here is not 100% correct, but the difference is small, for fairly collimated particles.
503 RecalculateMomentum(pjet,fPmissing);
504 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
506 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
508 ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
509 Double_t ptdiff=fPtJet-origPtJet;
510 ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
511 ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/origPtJet);
513 FillHistogramsRecoJetCorr(charm, jet, aodEvent);
515 }//end loop on candidates
517 //retrieve side band background candidates for Dstar
518 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
520 for(Int_t ib=0;ib<nsbcand;ib++){
521 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
522 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
523 if(!sbcand) continue;
525 fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
528 RecalculateMomentum(pjet,fPmissing);
529 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
531 SideBandBackground(sbcand,jet);
535 AliAODRecoDecayHF* charmbg = 0x0;
536 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
537 if(!charmbg) continue;
538 fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
542 RecalculateMomentum(pjet,fPmissing);
543 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
545 MCBackground(charmbg);
551 hNJetPerEv->Fill(cntjet);
552 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
557 //_______________________________________________________________________________
559 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
561 // The Terminate() function is the last function to be called during
562 // a query. It always runs on the client, it can be used to present
563 // the results graphically or save the results to file.
565 Info("Terminate"," terminate");
566 AliAnalysisTaskSE::Terminate();
568 fOutput = dynamic_cast<TList*> (GetOutputData(1));
570 printf("ERROR: fOutput not available\n");
575 //_______________________________________________________________________________
577 void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
579 Int_t abspdg=TMath::Abs(pdg);
581 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
582 // compute the Delta mass for the D*
583 if(fCandidateType==kDstartoKpipi){
585 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
589 fMinMass = mass-range;
590 fMaxMass = mass+range;
592 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
593 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
596 //_______________________________________________________________________________
598 void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
605 printf("Error! Lower limit larger than upper limit!\n");
606 fMinMass = uplimit - uplimit*0.2;
611 //_______________________________________________________________________________
613 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
615 AliInfo("Maximum number of bins allowed is 30!");
618 if(!width) return kFALSE;
619 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
623 //_______________________________________________________________________________
625 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
627 printf("Particle or jet do not exist!\n");
631 Bool_t okpp=part->PxPyPz(p);
632 Bool_t okpj=jet->PxPyPz(pj);
634 printf("Problems getting momenta\n");
638 RecalculateMomentum(pj, fPmissing);
639 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
640 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
644 //_______________________________________________________________________________
646 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
654 //_______________________________________________________________________________
656 Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
659 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
660 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
661 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
662 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
663 hstat->GetXaxis()->SetBinLabel(4,"N jets");
664 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
665 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
666 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
667 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
668 hstat->SetNdivisions(1);
671 const Int_t nbinsmass=200;
672 const Int_t nbinsptjet=500;
673 const Int_t nbinsptD=100;
674 const Int_t nbinsz=100;
675 const Int_t nbinsphi=200;
676 const Int_t nbinseta=100;
677 const Int_t nbinsContrib=100;
678 const Int_t nbinsA=100;
680 const Float_t ptjetlims[2]={0.,200.};
681 const Float_t ptDlims[2]={0.,50.};
682 const Float_t zlims[2]={0.,1.2};
683 const Float_t philims[2]={0.,6.3};
684 const Float_t etalims[2]={-1.5,1.5};
685 const Int_t nContriblims[2]={0,100};
686 const Float_t arealims[2]={0.,2};
688 // jet related fistograms
690 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
692 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
693 hPhiJetTrks->Sumw2();
694 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
695 hEtaJetTrks->Sumw2();
696 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
699 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
701 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
703 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
705 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
708 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
709 hdeltaRJetTracks->Sumw2();
711 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
714 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
718 fOutput->Add(hEjetTrks);
719 fOutput->Add(hPhiJetTrks);
720 fOutput->Add(hEtaJetTrks);
721 fOutput->Add(hPtJetTrks);
723 fOutput->Add(hPhiJet);
724 fOutput->Add(hEtaJet);
725 fOutput->Add(hPtJet);
726 fOutput->Add(hdeltaRJetTracks);
727 fOutput->Add(hNtrArr);
728 fOutput->Add(hNJetPerEv);
733 const Int_t nbinsSparse[nAxis]={nbinsphi,nbinseta, nbinsptjet, nbinsptjet,nbinsContrib, nbinsA};
734 const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],nContriblims[0],arealims[0]};
736 maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],nContriblims[1],arealims[1]};
737 THnSparseF *hsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
745 TH1I* hControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
746 hControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
747 hControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
748 hControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
749 hControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
750 hControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
751 hControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
752 hControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
753 hControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
755 hControlDInJ->SetNdivisions(1);
756 hControlDInJ->GetXaxis()->SetLabelSize(0.05);
757 fOutput->Add(hControlDInJ);
759 TH1F *hmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo)", 50, 0.,30.);
760 fOutput->Add(hmissingp);
761 TH1F *hDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction", 50, 0.,30.);
762 fOutput->Add(hDeltaPtJet);
763 TH1F *hRelDeltaPtJet=new TH1F("hRelDeltaPtJet", "Difference between the jet pt before and after correction/ original pt", 200, 0.,2.);
764 fOutput->Add(hRelDeltaPtJet);
766 TH1I* hIDddaugh=new TH1I("hIDddaugh", "ID of daughters", 2001,-1000,1000);
767 fOutput->Add(hIDddaugh);
768 TH1I* hIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet", 2001,-1000,1000);
769 fOutput->Add(hIDddaughOut);
770 TH1I* hIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks", 2001,-1000,1000);
771 fOutput->Add(hIDjetTracks);
773 if(fCandidateType==kDstartoKpipi)
776 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
777 hDiffSideBand->SetStats(kTRUE);
778 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
779 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
780 hDiffSideBand->Sumw2();
781 fOutput->Add(hDiffSideBand);
784 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
785 hPtPion->SetStats(kTRUE);
786 hPtPion->GetXaxis()->SetTitle("GeV/c");
787 hPtPion->GetYaxis()->SetTitle("Entries");
789 fOutput->Add(hPtPion);
792 // D related histograms
793 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
794 hNDPerEvNoJet->Sumw2();
795 fOutput->Add(hNDPerEvNoJet);
797 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
798 hptDPerEvNoJet->Sumw2();
799 fOutput->Add(hptDPerEvNoJet);
801 const Int_t nAxisD=4;
802 const Int_t nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
803 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
804 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
805 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
806 hsDstandalone->Sumw2();
808 fOutput->Add(hsDstandalone);
810 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]);
811 hPtJetWithD->Sumw2();
812 //for the MC this histogram is filled with the real background
813 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]);
814 hPtJetWithDsb->Sumw2();
816 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);
817 hNJetPerEvNoD->Sumw2();
819 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
820 hPtJetPerEvNoD->Sumw2();
822 fOutput->Add(hPtJetWithD);
823 fOutput->Add(hPtJetWithDsb);
824 fOutput->Add(hNJetPerEvNoD);
825 fOutput->Add(hPtJetPerEvNoD);
827 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
829 fOutput->Add(hDeltaRD);
831 //background (side bands for the Dstar and like sign for D0)
832 fJetRadius=GetJetContainer(0)->GetJetRadius();
833 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]);
834 hInvMassptD->SetStats(kTRUE);
835 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
836 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
837 hInvMassptD->Sumw2();
839 fOutput->Add(hInvMassptD);
842 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]);
843 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
844 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
845 hInvMassptDbg->Sumw2();
846 fOutput->Add(hInvMassptDbg);
849 Double_t pi=TMath::Pi(), philims2[2];
850 philims2[0]=-pi*1./2.;
851 philims2[1]=pi*3./2.;
853 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
854 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
855 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
856 THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. For SB or not and within the jet cone or not", nAxis, nbinsSparse, minSparse, maxSparse);
859 fOutput->Add(hsDphiz);
861 PostData(1, fOutput);
866 //_______________________________________________________________________________
868 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet, AliAODEvent* aodEvent){
870 Double_t ptD=candidate->Pt();
871 Double_t deltaR=DeltaR(candidate,jet);
872 Double_t phiD=candidate->Phi();
873 Double_t deltaphi = jet->Phi()-phiD;
874 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
875 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
876 Double_t z=Z(candidate,jet);
877 if(z>1) ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
879 TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
880 hDeltaRD->Fill(deltaR);
882 if(fCandidateType==kD0toKpi) {
883 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
885 FillHistogramsD0JetCorr(dzero,deltaphi,phiD,z,ptD,fPtJet, aodEvent);
889 if(fCandidateType==kDstartoKpipi) {
890 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
891 FillHistogramsDstarJetCorr(dstar,deltaphi,phiD,z,ptD,fPtJet);
894 } else{ //generated level
895 //AliInfo("Non reco");
896 FillHistogramsMCGenDJetCorr(deltaphi,phiD,z,ptD,fPtJet);
900 //_______________________________________________________________________________
902 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t phiD, Double_t z, Double_t ptD, Double_t ptj, AliAODEvent* aodEvent){
905 Double_t masses[2]={0.,0.};
906 Int_t pdgdaughtersD0[2]={211,321};//pi,K
907 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
909 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
910 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
912 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
913 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
914 Double_t point[8]={z,dPhi,ptj,ptD,masses[0],0, fIsDInJet ? 1 : 0,phiD};
915 Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
916 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
917 if(isselected==1 || isselected==3) {
919 if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
921 FillMassHistograms(masses[0], ptD);
922 hsDphiz->Fill(point,1.);
925 if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
927 FillMassHistograms(masses[1], ptD);
929 hsDphiz->Fill(point,1.);
934 //_______________________________________________________________________________
936 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t phiD, Double_t z, Double_t ptD, Double_t ptj){
937 //dPhi and z not used at the moment,but will be (re)added
939 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
940 Double_t deltamass= dstar->DeltaInvMass();
941 //Double_t massD0= dstar->InvMassD0();
944 TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
945 hPtPion->Fill(softpi->Pt());
947 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
948 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
949 Int_t isSB=0;//IsDzeroSideBand(dstar);
951 Double_t point[8]={z,dPhi,ptj,ptD,deltamass,isSB, fIsDInJet ? 1 : 0,phiD};
953 if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
955 FillMassHistograms(deltamass, ptD);
956 hsDphiz->Fill(point,1.);
960 //_______________________________________________________________________________
962 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t phiD,Double_t z,Double_t ptD,Double_t ptjet){
965 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
966 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
967 Double_t point[8]={z,dPhi,ptjet,ptD,pdgmass,0, fIsDInJet ? 1 : 0,phiD};
969 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
970 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
972 hsDphiz->Fill(point,1.);
974 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
979 //_______________________________________________________________________________
981 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
984 TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
985 hInvMassptD->Fill(mass,ptD);
989 //________________________________________________________________________________
991 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
993 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
994 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
995 if (fIsDInJet) jet->AddFlavourTag(tag);
1001 //_______________________________________________________________________________
1003 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1005 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1006 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1007 // has a width of ~5 sigmas. Two band needed for opening angle considerations
1008 TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1009 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1010 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1012 Double_t deltaM=candDstar->DeltaInvMass();
1013 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1014 Double_t z=Z(candDstar,jet);
1015 Double_t ptD=candDstar->Pt();
1017 Double_t dPhi=jet->Phi()-candDstar->Phi();
1019 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1020 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1022 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1023 Double_t point[7]={z,dPhi,fPtJet,ptD,deltaM,isSideBand, fIsDInJet ? 1 : 0};
1024 //fill the background histogram with the side bands or when looking at MC with the real background
1026 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1027 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1028 hsDphiz->Fill(point,1.);
1029 if(fIsDInJet){ // evaluate in the near side
1030 hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1035 //_______________________________________________________________________________
1037 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg){
1039 //need mass, deltaR, pt jet, ptD
1040 //two cases: D0 and Dstar
1042 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1043 if(!isselected) return;
1045 Double_t ptD=candbg->Pt();
1047 TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1048 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1051 if(fCandidateType==kDstartoKpipi){
1052 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1053 Double_t deltaM=dstarbg->DeltaInvMass();
1054 hInvMassptDbg->Fill(deltaM,ptD);
1055 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1058 if(fCandidateType==kD0toKpi){
1059 Double_t masses[2]={0.,0.};
1060 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1061 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1063 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1064 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1067 if(isselected==1 || isselected==3) {
1068 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1069 hInvMassptDbg->Fill(masses[0],ptD);
1072 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1073 hInvMassptDbg->Fill(masses[1],ptD);
1080 //_______________________________________________________________________________
1082 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1083 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1085 if(!p1 || !p2) return -1;
1086 Double_t phi1=p1->Phi(),eta1=p1->Eta();
1087 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1089 Double_t dPhi=phi1-phi2;
1090 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1091 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1093 Double_t dEta=eta1-eta2;
1094 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1099 //_______________________________________________________________________________
1101 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1103 Double_t ptD=candDstar->Pt();
1104 Int_t bin = fCuts->PtBin(ptD);
1106 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1107 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1108 return -1; // out of bounds
1111 Double_t invM=candDstar->InvMassD0();
1112 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1114 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1115 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1117 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1122 //_______________________________________________________________________________
1124 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1125 //returns true/false and the array of particles not found among jet constituents
1127 TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1128 TH1I* hIDddaugh =(TH1I*)fOutput->FindObject("hIDddaugh");
1129 TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1130 TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1132 Int_t daughtersID[3];
1134 Int_t dmatchedID[3]={0,0,0};
1135 Int_t countmatches=0;
1136 if (fCandidateType==kDstartoKpipi){
1137 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1138 AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1139 daughtersID[0]=dzero->GetProngID(0);
1140 daughtersID[1]=dzero->GetProngID(1);
1141 daughtersID[2]=dstar->GetBachelor()->GetID();
1144 dtrks=new AliAODTrack*[3];
1145 dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1146 dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1147 dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1151 if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID()) hControlDInJ->Fill(6);
1153 hIDddaugh->Fill(daughtersID[0]);
1154 hIDddaugh->Fill(daughtersID[1]);
1155 hIDddaugh->Fill(daughtersID[2]);
1158 //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1161 if (fCandidateType==kD0toKpi){
1162 daughtersID[0]=charm->GetProngID(0);
1163 daughtersID[1]=charm->GetProngID(1);
1166 hIDddaugh->Fill(daughtersID[0]);
1167 hIDddaugh->Fill(daughtersID[1]);
1169 dtrks=new AliAODTrack*[2];
1170 dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1171 dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1175 const Int_t cndaugh=ndaugh;
1176 daughOutOfJetID=new Int_t[cndaugh];
1178 Int_t nchtrk=thejet->GetNumberOfTracks();
1179 for(Int_t itk=0;itk<nchtrk;itk++){
1180 AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1181 Int_t idtkjet=tkjet->GetID();
1182 if(fillH) hIDjetTracks->Fill(idtkjet);
1183 for(Int_t id=0;id<ndaugh;id++){
1184 if(idtkjet==daughtersID[id]) {
1185 dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1190 if(countmatches==ndaugh) break;
1193 //reverse: include in the array the ID of daughters not matching
1195 for(Int_t id=0;id<ndaugh;id++){
1196 if(dmatchedID[id]==0) {
1197 daughOutOfJetID[id]=daughtersID[id];
1198 if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1200 else daughOutOfJetID[id]=0;
1203 if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1204 if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1205 if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1207 if(ndaugh!=countmatches){
1214 //_______________________________________________________________________________
1216 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1218 //check the conditions type:
1219 //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1220 //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1221 //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1223 TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1225 Bool_t testDeltaR=kFALSE;
1226 if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1228 Int_t* daughOutOfJet;
1229 AliAODTrack** charmDaugh;
1230 Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1231 if(!testDaugh && testDeltaR && fTypeDInJet==2){
1234 if(fCandidateType==kD0toKpi) ndaugh=2;
1240 for(Int_t id=0;id<ndaugh;id++){
1241 if(daughOutOfJet[id]!=0){ //non-matched daughter
1242 //get track and its momentum
1245 hControlDInJ->Fill(3);
1246 if(id<2) hControlDInJ->Fill(4);
1247 if(id==2)hControlDInJ->Fill(5);
1249 fPmissing[0]+=charmDaugh[id]->Px();
1250 fPmissing[1]+=charmDaugh[id]->Py();
1251 fPmissing[2]+=charmDaugh[id]->Pz();
1256 //now the D in within the jet
1261 delete[] charmDaugh;
1264 switch(fTypeDInJet){
1268 result=testDeltaR && testDaugh;
1270 result=testDeltaR && testDaugh;