]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
dont create any files in place to avoid problems with eos
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskFlavourJetCorrelations.cxx
CommitLineData
931f8b00 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// A. Grelli (Utrecht University) a.grelli@uu.nl
25// X. Zhang (LBNL) XMZhang@lbl.gov
26//-----------------------------------------------------------------------
27
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include "TROOT.h"
31#include <TH3F.h>
32#include <THnSparse.h>
33
34#include "AliAnalysisTaskFlavourJetCorrelations.h"
35#include "AliAODMCHeader.h"
36#include "AliAODHandler.h"
37#include "AliAnalysisManager.h"
38#include "AliLog.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"
49
50ClassImp(AliAnalysisTaskFlavourJetCorrelations)
51
52
53//_______________________________________________________________________________
54
55AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
0455d618 56AliAnalysisTaskEmcalJet("",kTRUE),
931f8b00 57fUseMCInfo(kTRUE),
58fUseReco(kTRUE),
59fCandidateType(),
60fPDGmother(),
61fNProngs(),
62fPDGdaughters(),
63fBranchName(),
931f8b00 64fCuts(0),
65fMinMass(),
66fMaxMass(),
67fJetArrName(0),
68fCandArrName(0),
69fLeadingJetOnly(kFALSE),
bbb94467 70fJetRadius(0),
71fCandidateArray(0),
0455d618 72fSideBandArray(0),
b2705b43 73fJetOnlyMode(0),
74fPmissing(),
75fPtJet(0),
76fIsDInJet(0),
77fTypeDInJet(2),
78fTrackArr(0)
931f8b00 79{
80 //
81 // Default ctor
82 //
83 //SetMakeGeneralHistograms(kTRUE);
84
85}
86
87//_______________________________________________________________________________
88
a8b2b672 89AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, Bool_t jetOnly) :
0455d618 90AliAnalysisTaskEmcalJet(name,kTRUE),
931f8b00 91fUseMCInfo(kTRUE),
92fUseReco(kTRUE),
93fCandidateType(),
94fPDGmother(),
95fNProngs(),
96fPDGdaughters(),
97fBranchName(),
931f8b00 98fCuts(0),
99fMinMass(),
100fMaxMass(),
101fJetArrName(0),
102fCandArrName(0),
103fLeadingJetOnly(kFALSE),
bbb94467 104fJetRadius(0),
105fCandidateArray(0),
0455d618 106fSideBandArray(0),
b2705b43 107fJetOnlyMode(jetOnly),
108fPmissing(),
109fPtJet(0),
110fIsDInJet(0),
111fTypeDInJet(2),
112fTrackArr(0)
931f8b00 113{
114 //
115 // Constructor. Initialization of Inputs and Outputs
116 //
117
118 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
119 fCuts=cuts;
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};
123
124 switch(fCandidateType){
125 case 0:
126 fPDGmother=421;
127 fNProngs=2;
128 fPDGdaughters[0]=211;//pi
129 fPDGdaughters[1]=321;//K
130 fPDGdaughters[2]=0; //empty
131 fPDGdaughters[3]=0; //empty
132 fBranchName="D0toKpi";
133 fCandArrName="D0";
134 break;
135 case 1:
136 fPDGmother=413;
137 fNProngs=3;
138 fPDGdaughters[1]=211;//pi soft
139 fPDGdaughters[0]=421;//D0
140 fPDGdaughters[2]=211;//pi fromD0
141 fPDGdaughters[3]=321; // K from D0
142 fBranchName="Dstar";
143 fCandArrName="Dstar";
144
145 if(nptbins<=13){
146 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
147 } else {
148 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
149 }
150 break;
151 default:
152 printf("%d not accepted!!\n",fCandidateType);
153 break;
154 }
155
156 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
157 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
158
0455d618 159 if(fJetOnlyMode){
160 DefineOutput(1,TList::Class()); //histos with jet info
161 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
162 }
163 else{
164 DefineInput(1, TClonesArray::Class());
165 DefineInput(2, TClonesArray::Class());
166
167 DefineOutput(1,TList::Class()); // histos
168 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
169 }
931f8b00 170}
171
172//_______________________________________________________________________________
173
174AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
175 //
176 // destructor
177 //
178
179 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
180
931f8b00 181 delete fCuts;
182
183}
184
185//_______________________________________________________________________________
186
187void AliAnalysisTaskFlavourJetCorrelations::Init(){
188 //
189 // Initialization
190 //
191
192 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
0455d618 193
931f8b00 194 switch(fCandidateType){
195 case 0:
196 {
197 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
198 copyfCuts->SetName("AnalysisCutsDzero");
199 // Post the data
200 PostData(2,copyfCuts);
201 }
202 break;
203 case 1:
204 {
205 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
206 copyfCuts->SetName("AnalysisCutsDStar");
207 // Post the cuts
208 PostData(2,copyfCuts);
209 }
210 break;
211 default:
212 return;
213 }
214
215 return;
216}
217
218//_______________________________________________________________________________
219
220void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
221 // output
222 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
76bf81f2 223 AliAnalysisTaskEmcal::UserCreateOutputObjects();
931f8b00 224 // define histograms
76bf81f2 225 // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
931f8b00 226 DefineHistoForAnalysis();
76bf81f2 227 PostData(1,fOutput);
931f8b00 228
229 return;
230}
231
232//_______________________________________________________________________________
233
b5d0ba13 234Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
931f8b00 235{
76bf81f2 236 // user exec from AliAnalysisTaskEmcal is used
237
931f8b00 238 // Load the event
239 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
240
241 TClonesArray *arrayDStartoD0pi=0;
931f8b00 242
243 if (!aodEvent && AODEvent() && IsStandardAOD()) {
244
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());
248
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());
257 }
258 } else {
259 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
260 }
261
262 if (!arrayDStartoD0pi) {
263 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
264 // return;
265 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
266
267 TClonesArray* mcArray = 0x0;
268 if (fUseMCInfo) {
269 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
270 if (!mcArray) {
271 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
b5d0ba13 272 return kFALSE;
931f8b00 273 }
274 }
275
276 //retrieve jets
b2705b43 277 fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
931f8b00 278 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
279 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
280 fJetRadius=GetJetContainer(0)->GetJetRadius();
281
b2705b43 282 if(!fTrackArr){
931f8b00 283 AliInfo(Form("Could not find the track array\n"));
b5d0ba13 284 return kFALSE;
931f8b00 285 }
286
bbb94467 287
c3c1ab31 288 fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
b5d0ba13 289 if (!fCandidateArray) return kFALSE;
c3c1ab31 290 if (fCandidateType==1) {
291 fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
b5d0ba13 292 if (!fSideBandArray) return kFALSE;
c3c1ab31 293 }
294 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
931f8b00 295
296 //Histograms
76bf81f2 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");
0455d618 314 THnSparseF* hsJet=(THnSparseF*)fOutput->FindObject("hsJet");
315
09fce7b3 316 TH1F* hztracksinjet=(TH1F*)fOutput->FindObject("hztracksinjet");
317 TH1F* hDiffPtTrPtJzNok=(TH1F*)fOutput->FindObject("hDiffPtTrPtJzNok");
318 TH1F* hDiffPtTrPtJzok=(TH1F*)fOutput->FindObject("hDiffPtTrPtJzok");
319 TH1F* hDiffPzTrPtJzok=(TH1F*)fOutput->FindObject("hDiffPzTrPtJzok");
320 TH1F* hDiffPzTrPtJzNok=(TH1F*)fOutput->FindObject("hDiffPzTrPtJzNok");
321 TH1I* hNtrkjzNok=(TH1I*)fOutput->FindObject("hNtrkjzNok");
322 TH1F* hztracksinjetT=(TH1F*)fOutput->FindObject("hztracksinjetT");
323
324
931f8b00 325 hstat->Fill(0);
326
327 // fix for temporary bug in ESDfilter
328 // the AODs with null vertex pointer didn't pass the PhysSel
b5d0ba13 329 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
931f8b00 330
331 //Event selection
332 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
333 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
b5d0ba13 334 if(!iseventselected) return kFALSE;
931f8b00 335
336 hstat->Fill(1);
931f8b00 337
338 //retrieve charm candidates selected
0455d618 339 Int_t candidates = 0;
d9126a4c 340 Int_t njets=GetJetContainer()->GetNJets();
0455d618 341
342 if(!fJetOnlyMode) {
343 candidates = fCandidateArray->GetEntriesFast();
344
345 //trigger on jets
d9126a4c
CB		 346 if(njets == 0) {
347 hstat->Fill(6, candidates);
348 hNDPerEvNoJet->Fill(candidates);
349 for(Int_t iD=0;iD<candidates;iD++){
bbb94467 350 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
d9126a4c
CB		 351 if(!cand) continue;
352 hptDPerEvNoJet->Fill(cand->Pt());
353
354 }
b5d0ba13 355 return kFALSE;
d9126a4c
CB		 356
357 }
1353fa73 358
359 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
360
361 for(Int_t ic = 0; ic < candidates; ic++) {
362
363 // D* candidates
bbb94467 364 AliAODRecoDecayHF* charm=0x0;
365 AliAODRecoCascadeHF* dstar=0x0;
366
367
368 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
1353fa73 369 if(!charm) continue;
bbb94467 370
371 if(fCandidateType==kDstartoKpipi){
372 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
373 if(!dstar) continue;
374 }
375
1353fa73 376 hstat->Fill(2);
377
378 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
379
380 if(fCandidateType==kDstartoKpipi) {
8577f7bc 381 if(fUseReco){
382 Double_t deltamass= dstar->DeltaInvMass();
383 candsparse[3]=deltamass;
384 } else candsparse[3] = 0.145; //for generated
1353fa73 385 hnspDstandalone->Fill(candsparse);
386 }
387 if(fCandidateType==kD0toKpi){
8577f7bc 388 if(fUseReco){
389 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
390 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
391 //not a further selection,just to get the value of isselected for the D0
392 Double_t masses[2];
393 Int_t pdgdaughtersD0[2]={211,321};//pi,K
394 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
395
396 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
397 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
398 if(isselected==1 || isselected==3) {
399 candsparse[3]=masses[0];
400 hnspDstandalone->Fill(candsparse);
401 }
402 if(isselected>=2){
403 candsparse[3]=masses[1];
404 hnspDstandalone->Fill(candsparse);
405
406 }
407 } else { //generated
408 Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
409 candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1353fa73 410 hnspDstandalone->Fill(candsparse);
1353fa73 411 }
412 }
413 }
0455d618 414 }
415
931f8b00 416 // we start with jets
417 Double_t ejet = 0;
418 Double_t phiJet = 0;
419 Double_t etaJet = 0;
420 Double_t ptjet = 0;
421 Double_t leadingJet =0;
0455d618 422 Double_t pointJ[6];
931f8b00 423
b2705b43 424 Int_t ntrarr=fTrackArr->GetEntriesFast();
931f8b00 425 hNtrArr->Fill(ntrarr);
426
427 for(Int_t i=0;i<ntrarr;i++){
b2705b43 428 AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
931f8b00 429 //reusing histograms
430 hPtJetTrks->Fill(jtrack->Pt());
431 hPhiJetTrks->Fill(jtrack->Phi());
432 hEtaJetTrks->Fill(jtrack->Eta());
433 hEjetTrks->Fill(jtrack->E());
434 }
435
436
437 //option to use only the leading jet
438 if(fLeadingJetOnly){
439 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
440 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
441 ptjet = jetL->Pt();
442 if(ptjet>leadingJet ) leadingJet = ptjet;
443
444 }
445 }
446
447 Int_t cntjet=0;
448 //loop over jets and consider only the leading jet in the event
09fce7b3 449 for (Int_t iJets = 0; iJets<njets; iJets++) {
450 fPmissing[0]=0;
451 fPmissing[1]=0;
452 fPmissing[2]=0;
453
931f8b00 454 //Printf("Jet N %d",iJets);
455 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
456 //Printf("Corr task Accept Jet");
457 if(!AcceptJet(jet)) {
458 hstat->Fill(5);
459 continue;
460 }
461
462 //jets variables
463 ejet = jet->E();
464 phiJet = jet->Phi();
465 etaJet = jet->Eta();
b2705b43 466 fPtJet = jet->Pt();
467 Double_t origPtJet=fPtJet;
09fce7b3 468
0455d618 469 pointJ[0] = phiJet;
470 pointJ[1] = etaJet;
471 pointJ[2] = ptjet;
472 pointJ[3] = ejet;
473 pointJ[4] = jet->GetNumberOfConstituents();
474 pointJ[5] = jet->Area();
931f8b00 475
476 // choose the leading jet
477 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
478 //used for normalization
479 hstat->Fill(3);
480 cntjet++;
481 // fill energy, eta and phi of the jet
482 hEjet ->Fill(ejet);
483 hPhiJet ->Fill(phiJet);
484 hEtaJet ->Fill(etaJet);
b2705b43 485 hPtJet ->Fill(fPtJet);
0455d618 486 if(fJetOnlyMode) hsJet->Fill(pointJ,1);
931f8b00 487 //loop on jet particles
09fce7b3 488 Int_t ntrjet= jet->GetNumberOfTracks();
489 Double_t sumPtTracks=0,sumPzTracks=0;
490 Int_t zg1jtrk=0;
931f8b00 491 for(Int_t itrk=0;itrk<ntrjet;itrk++){
492
b2705b43 493 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
931f8b00 494 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
09fce7b3 495 Double_t ztrackj=Z(jetTrk,jet);
496 hztracksinjet->Fill(ztrackj);
497 sumPtTracks+=jetTrk->Pt();
498 sumPzTracks+=jetTrk->Pz();
499 if(ztrackj>1){
500 zg1jtrk++;
501 }
502
503 Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
504 hztracksinjetT->Fill(ztrackjTr);
505
931f8b00 506 }//end loop on jet tracks
507
09fce7b3 508 hNtrkjzNok->Fill(zg1jtrk);
509 Double_t diffpt=origPtJet-sumPtTracks;
510 Double_t diffpz=jet->Pz()-sumPzTracks;
511 if(zg1jtrk>0){
512 hDiffPtTrPtJzNok->Fill(diffpt);
513 hDiffPzTrPtJzNok->Fill(diffpz);
514
515 }else{
516 hDiffPtTrPtJzok->Fill(diffpt);
517 hDiffPzTrPtJzok->Fill(diffpz);
518 }
519
d9126a4c
CB		 520 if(candidates==0){
521 hstat->Fill(7);
b2705b43 522 hPtJetPerEvNoD->Fill(fPtJet);
d9126a4c 523 }
0455d618 524 if(!fJetOnlyMode) {
b2705b43 525 //Printf("N candidates %d ", candidates);
526 for(Int_t ic = 0; ic < candidates; ic++) {
527
528 // D* candidates
529 AliVParticle* charm=0x0;
530 charm=(AliVParticle*)fCandidateArray->At(ic);
531 if(!charm) continue;
532 AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
533 fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
534 if (fIsDInJet) FlagFlavour(jet);
535 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
536
537 //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.
538
539 Double_t pjet[3];
540 jet->PxPyPz(pjet);
541 RecalculateMomentum(pjet,fPmissing);
542 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
543
09fce7b3 544
545 //debugging histograms
b2705b43 546 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
09fce7b3 547 for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
b2705b43 548
549 ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
550 Double_t ptdiff=fPtJet-origPtJet;
551 ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
552 ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/origPtJet);
553
554 FillHistogramsRecoJetCorr(charm, jet, aodEvent);
555
556 }//end loop on candidates
931f8b00 557
b2705b43 558 //retrieve side band background candidates for Dstar
559 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
931f8b00 560
b2705b43 561 for(Int_t ib=0;ib<nsbcand;ib++){
562 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
563 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
564 if(!sbcand) continue;
565
566 fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
567 Double_t pjet[3];
568 jet->PxPyPz(pjet);
569 RecalculateMomentum(pjet,fPmissing);
570 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
571
572 SideBandBackground(sbcand,jet);
573
574 }
575 if(fUseMCInfo){
576 AliAODRecoDecayHF* charmbg = 0x0;
577 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
578 if(!charmbg) continue;
579 fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
580
581 Double_t pjet[3];
582 jet->PxPyPz(pjet);
583 RecalculateMomentum(pjet,fPmissing);
584 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
585
586 MCBackground(charmbg);
587 }
931f8b00 588 }
589 }
590 } // end of jet loop
591
592 hNJetPerEv->Fill(cntjet);
d9126a4c 593 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
76bf81f2 594 PostData(1,fOutput);
b5d0ba13 595 return kTRUE;
931f8b00 596}
597
598//_______________________________________________________________________________
599
600void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
601{
602 // The Terminate() function is the last function to be called during
603 // a query. It always runs on the client, it can be used to present
604 // the results graphically or save the results to file.
605
606 Info("Terminate"," terminate");
607 AliAnalysisTaskSE::Terminate();
608
76bf81f2 609 fOutput = dynamic_cast<TList*> (GetOutputData(1));
610 if (!fOutput) {
611 printf("ERROR: fOutput not available\n");
931f8b00 612 return;
613 }
614}
615
616//_______________________________________________________________________________
617
618void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
619 Float_t mass=0;
620 Int_t abspdg=TMath::Abs(pdg);
621
622 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
623 // compute the Delta mass for the D*
624 if(fCandidateType==kDstartoKpipi){
625 Float_t mass1=0;
626 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
627 mass = mass-mass1;
628 }
629
630 fMinMass = mass-range;
631 fMaxMass = mass+range;
632
633 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
634 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
635}
636
637//_______________________________________________________________________________
638
639void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
640 if(uplimit>lowlimit)
641 {
642 fMinMass = lowlimit;
643 fMaxMass = uplimit;
644 }
645 else{
646 printf("Error! Lower limit larger than upper limit!\n");
647 fMinMass = uplimit - uplimit*0.2;
648 fMaxMass = uplimit;
649 }
650}
651
652//_______________________________________________________________________________
653
654Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
655 if(nptbins>30) {
656 AliInfo("Maximum number of bins allowed is 30!");
657 return kFALSE;
658 }
659 if(!width) return kFALSE;
660 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
661 return kTRUE;
662}
663
664//_______________________________________________________________________________
665
09fce7b3 666Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
931f8b00 667 if(!part || !jet){
668 printf("Particle or jet do not exist!\n");
669 return -99;
670 }
671 Double_t p[3],pj[3];
672 Bool_t okpp=part->PxPyPz(p);
673 Bool_t okpj=jet->PxPyPz(pj);
674 if(!okpp || !okpj){
675 printf("Problems getting momenta\n");
676 return -999;
677 }
b2705b43 678
679 RecalculateMomentum(pj, fPmissing);
09fce7b3 680 if(transverse) return ZT(p,pj);
681 else
682 return Z(p,pj);
683}
684
685//_______________________________________________________________________________
686Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
687
b2705b43 688 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
689 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
931f8b00 690 return z;
691}
692
09fce7b3 693
694//_______________________________________________________________________________
695Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
696
697 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
698 Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
699 return z;
700}
701
931f8b00 702//_______________________________________________________________________________
703
b2705b43 704void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
705
706 pj[0]+=pmissing[0];
707 pj[1]+=pmissing[1];
708 pj[2]+=pmissing[2];
709
710}
711
712//_______________________________________________________________________________
713
931f8b00 714Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
715
716 // Statistics
d9126a4c 717 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
931f8b00 718 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
719 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
d9126a4c 720 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
931f8b00 721 hstat->GetXaxis()->SetBinLabel(4,"N jets");
722 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
723 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
d9126a4c
CB		 724 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
725 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
931f8b00 726 hstat->SetNdivisions(1);
76bf81f2 727 fOutput->Add(hstat);
931f8b00 728
729 const Int_t nbinsmass=200;
730 const Int_t nbinsptjet=500;
731 const Int_t nbinsptD=100;
732 const Int_t nbinsz=100;
733 const Int_t nbinsphi=200;
1353fa73 734 const Int_t nbinseta=100;
0455d618 735 const Int_t nbinsContrib=100;
736 const Int_t nbinsA=100;
737
931f8b00 738 const Float_t ptjetlims[2]={0.,200.};
739 const Float_t ptDlims[2]={0.,50.};
740 const Float_t zlims[2]={0.,1.2};
741 const Float_t philims[2]={0.,6.3};
1353fa73 742 const Float_t etalims[2]={-1.5,1.5};
0455d618 743 const Int_t nContriblims[2]={0,100};
744 const Float_t arealims[2]={0.,2};
1353fa73 745
931f8b00 746 // jet related fistograms
747
748 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
749 hEjetTrks->Sumw2();
750 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
751 hPhiJetTrks->Sumw2();
1353fa73 752 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 753 hEtaJetTrks->Sumw2();
754 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
755 hPtJetTrks->Sumw2();
756
757 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
758 hEjet->Sumw2();
759 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
760 hPhiJet->Sumw2();
1353fa73 761 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 762 hEtaJet->Sumw2();
763 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
764 hPtJet->Sumw2();
765
931f8b00 766 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
767 hdeltaRJetTracks->Sumw2();
768
769 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
770 hNtrArr->Sumw2();
d9126a4c 771
931f8b00 772 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
773 hNJetPerEv->Sumw2();
774
d9126a4c 775
76bf81f2 776 fOutput->Add(hEjetTrks);
777 fOutput->Add(hPhiJetTrks);
778 fOutput->Add(hEtaJetTrks);
779 fOutput->Add(hPtJetTrks);
780 fOutput->Add(hEjet);
781 fOutput->Add(hPhiJet);
782 fOutput->Add(hEtaJet);
783 fOutput->Add(hPtJet);
76bf81f2 784 fOutput->Add(hdeltaRJetTracks);
785 fOutput->Add(hNtrArr);
786 fOutput->Add(hNJetPerEv);
931f8b00 787
0455d618 788 if(fJetOnlyMode){
789 //thnsparse for jets
790 const Int_t nAxis=6;
791 const Int_t nbinsSparse[nAxis]={nbinsphi,nbinseta, nbinsptjet, nbinsptjet,nbinsContrib, nbinsA};
2942f542 792 const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
0455d618 793 const Double_t
2942f542 794 maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
0455d618 795 THnSparseF *hsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
796 hsJet->Sumw2();
931f8b00 797
0455d618 798 fOutput->Add(hsJet);
799
931f8b00 800 }
931f8b00 801
0455d618 802 if(!fJetOnlyMode){
09fce7b3 803
804 //debugging histograms
805 TH1I* hControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
806 hControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
807 hControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
808 hControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
809 hControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
810 hControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
811 hControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
812 hControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
813 hControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
814
815 hControlDInJ->SetNdivisions(1);
816 hControlDInJ->GetXaxis()->SetLabelSize(0.05);
817 fOutput->Add(hControlDInJ);
818
819 TH1F *hmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
820 fOutput->Add(hmissingp);
821
822 for(Int_t k=0;k<3;k++) {
823 TH1F *hMissPi=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
824 fOutput->Add(hMissPi);
825 }
826 TH1F *hDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
827
828 fOutput->Add(hDeltaPtJet);
829 TH1F *hRelDeltaPtJet=new TH1F("hRelDeltaPtJet", "Difference between the jet pt before and after correction/ original pt;(p_{T}^{jet,recalc}-p_{T}^{jet,orig})/p_{T}^{jet,orig}", 200, 0.,1.);
830 fOutput->Add(hRelDeltaPtJet);
831
832 TH1F* hzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
833 fOutput->Add(hzDinjet);
834 //frag func of tracks in the jet
835 TH1F* hztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
836 fOutput->Add(hztracksinjet);
837
838 //check jet and tracks
839 TH1F* hDiffPtTrPtJzok=new TH1F("hDiffPtTrPtJzok","Sum p_{T}^{trks} - p_{T}^{jet}, for z<1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
840 fOutput->Add(hDiffPtTrPtJzok);
841 TH1F* hDiffPtTrPtJzNok=new TH1F("hDiffPtTrPtJzNok","Sum p_{T}^{trks} - p_{T}^{jet}, for z>1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
842 fOutput->Add(hDiffPtTrPtJzNok);
843 TH1F* hDiffPzTrPtJzok=new TH1F("hDiffPzTrPtJzok","Sum p_{z}^{trks} - p_{z}^{jet}, for z<1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
844 fOutput->Add(hDiffPzTrPtJzok);
845 TH1F* hDiffPzTrPtJzNok=new TH1F("hDiffPzTrPtJzNok","Sum p_{z}^{trks} - p_{z}^{jet}, for z>1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
846 fOutput->Add(hDiffPzTrPtJzNok);
847 TH1I* hNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
848 fOutput->Add(hNtrkjzNok);
849
850 //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
851 TH1F* hzDT=new TH1F("hzDT", "Z of D in jet in transverse components;(p_{T}^{D} dot p_{T}^{jet})/p_{T}^{jet}^{2} ",nbinsz,zlims[0],zlims[1]);
852 fOutput->Add(hzDT);
853 TH1F* hztracksinjetT=new TH1F("hztracksinjetT", "Z of jet tracks in transverse components;(p_{T}^{trks} dot p_{T}^{jet})/p_{T}^{jet}^{2}",nbinsz,zlims[0],zlims[1]);
854 fOutput->Add(hztracksinjetT);
855
856 TH1I* hIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
857 fOutput->Add(hIDddaugh);
858 TH1I* hIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
859 fOutput->Add(hIDddaughOut);
860 TH1I* hIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
861 fOutput->Add(hIDjetTracks);
862
863 TH1F* hDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
864 fOutput->Add(hDRdaughOut);
865
866
0455d618 867 if(fCandidateType==kDstartoKpipi)
868 {
869
870 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
871 hDiffSideBand->SetStats(kTRUE);
872 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
873 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
874 hDiffSideBand->Sumw2();
875 fOutput->Add(hDiffSideBand);
876
877
878 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
879 hPtPion->SetStats(kTRUE);
880 hPtPion->GetXaxis()->SetTitle("GeV/c");
881 hPtPion->GetYaxis()->SetTitle("Entries");
882 hPtPion->Sumw2();
883 fOutput->Add(hPtPion);
884
885 }
886 // D related histograms
887 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
888 hNDPerEvNoJet->Sumw2();
889 fOutput->Add(hNDPerEvNoJet);
890
891 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
892 hptDPerEvNoJet->Sumw2();
893 fOutput->Add(hptDPerEvNoJet);
894
895 const Int_t nAxisD=4;
896 const Int_t nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
897 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
898 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
899 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
900 hsDstandalone->Sumw2();
901
902 fOutput->Add(hsDstandalone);
903
904 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]);
905 hPtJetWithD->Sumw2();
906 //for the MC this histogram is filled with the real background
907 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]);
908 hPtJetWithDsb->Sumw2();
909
910 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);
911 hNJetPerEvNoD->Sumw2();
912
913 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
914 hPtJetPerEvNoD->Sumw2();
915
916 fOutput->Add(hPtJetWithD);
917 fOutput->Add(hPtJetWithDsb);
918 fOutput->Add(hNJetPerEvNoD);
919 fOutput->Add(hPtJetPerEvNoD);
920
921 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
922 hDeltaRD->Sumw2();
923 fOutput->Add(hDeltaRD);
924
925 //background (side bands for the Dstar and like sign for D0)
926 fJetRadius=GetJetContainer(0)->GetJetRadius();
927 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]);
928 hInvMassptD->SetStats(kTRUE);
929 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
930 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
931 hInvMassptD->Sumw2();
932
933 fOutput->Add(hInvMassptD);
934
935 if(fUseMCInfo){
936 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]);
937 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
938 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
939 hInvMassptDbg->Sumw2();
940 fOutput->Add(hInvMassptDbg);
941
942 }
943 Double_t pi=TMath::Pi(), philims2[2];
944 philims2[0]=-pi*1./2.;
945 philims2[1]=pi*3./2.;
946 const Int_t nAxis=7;
947 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
948 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
949 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
950 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);
951 hsDphiz->Sumw2();
952
953 fOutput->Add(hsDphiz);
954 }
76bf81f2 955 PostData(1, fOutput);
931f8b00 956
957 return kTRUE;
958}
959
960//_______________________________________________________________________________
961
8577f7bc 962void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet, AliAODEvent* aodEvent){
931f8b00 963
964 Double_t ptD=candidate->Pt();
931f8b00 965 Double_t deltaR=DeltaR(candidate,jet);
b2705b43 966 Double_t phiD=candidate->Phi();
967 Double_t deltaphi = jet->Phi()-phiD;
968 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
969 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
931f8b00 970 Double_t z=Z(candidate,jet);
09fce7b3 971 /*
972 if(z>1) {
973 ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
974 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
975
976 for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
977
978 ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
979 Double_t ptdiff=fPtJet-jet->Pt();
980 ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
981 ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
982
983
984 }
985 */
986 if(fIsDInJet)((TH1F*)fOutput->FindObject("hzDT"))->Fill(Z(candidate,jet,kTRUE));
931f8b00 987
76bf81f2 988 TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
931f8b00 989 hDeltaRD->Fill(deltaR);
990 if(fUseReco){
991 if(fCandidateType==kD0toKpi) {
992 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
8577f7bc 993
b2705b43 994 FillHistogramsD0JetCorr(dzero,deltaphi,phiD,z,ptD,fPtJet, aodEvent);
931f8b00 995
996 }
997
998 if(fCandidateType==kDstartoKpipi) {
999 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
b2705b43 1000 FillHistogramsDstarJetCorr(dstar,deltaphi,phiD,z,ptD,fPtJet);
931f8b00 1001
1002 }
1003 } else{ //generated level
1004 //AliInfo("Non reco");
b2705b43 1005 FillHistogramsMCGenDJetCorr(deltaphi,phiD,z,ptD,fPtJet);
931f8b00 1006 }
1007}
1008
1009//_______________________________________________________________________________
1010
b2705b43 1011void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t phiD, Double_t z, Double_t ptD, Double_t ptj, AliAODEvent* aodEvent){
931f8b00 1012
931f8b00 1013
1014 Double_t masses[2]={0.,0.};
1015 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1016 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1017
1018 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1019 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1020
76bf81f2 1021 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1022 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
2942f542 1023 Double_t point[8]={z,dPhi,ptj,ptD,masses[0],0, static_cast<Double_t>(fIsDInJet ? 1 : 0),phiD};
8577f7bc 1024 Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
931f8b00 1025 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1026 if(isselected==1 || isselected==3) {
1027
b2705b43 1028 if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
931f8b00 1029
b2705b43 1030 FillMassHistograms(masses[0], ptD);
931f8b00 1031 hsDphiz->Fill(point,1.);
1032 }
1033 if(isselected>=2) {
b2705b43 1034 if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
931f8b00 1035
b2705b43 1036 FillMassHistograms(masses[1], ptD);
931f8b00 1037 point[4]=masses[1];
1038 hsDphiz->Fill(point,1.);
1039 }
1040
1041}
1042
1043//_______________________________________________________________________________
1044
b2705b43 1045void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t phiD, Double_t z, Double_t ptD, Double_t ptj){
931f8b00 1046 //dPhi and z not used at the moment,but will be (re)added
1047
1048 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1049 Double_t deltamass= dstar->DeltaInvMass();
1050 //Double_t massD0= dstar->InvMassD0();
1051
1052
76bf81f2 1053 TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
931f8b00 1054 hPtPion->Fill(softpi->Pt());
1055
76bf81f2 1056 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1057 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
ed5a88ff 1058 Int_t isSB=0;//IsDzeroSideBand(dstar);
931f8b00 1059
2942f542 1060 Double_t point[8]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),phiD};
931f8b00 1061
b2705b43 1062 if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
931f8b00 1063
b2705b43 1064 FillMassHistograms(deltamass, ptD);
931f8b00 1065 hsDphiz->Fill(point,1.);
1066
1067}
1068
1069//_______________________________________________________________________________
1070
b2705b43 1071void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t phiD,Double_t z,Double_t ptD,Double_t ptjet){
931f8b00 1072
1073 Double_t pdgmass=0;
76bf81f2 1074 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1075 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
2942f542 1076 Double_t point[8]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),phiD};
931f8b00 1077
1078 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1079 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1080 point[4]=pdgmass;
8577f7bc 1081 hsDphiz->Fill(point,1.);
b2705b43 1082 if(fIsDInJet) {
931f8b00 1083 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
931f8b00 1084 }
1085
1086}
1087
1088//_______________________________________________________________________________
1089
b2705b43 1090void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
931f8b00 1091
b2705b43 1092 if(fIsDInJet) {
76bf81f2 1093 TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
931f8b00 1094 hInvMassptD->Fill(mass,ptD);
1095 }
1096}
1097
1353fa73 1098//________________________________________________________________________________
1099
b2705b43 1100void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1101
931f8b00 1102 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1103 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
b2705b43 1104 if (fIsDInJet) jet->AddFlavourTag(tag);
931f8b00 1105
1106 return;
1107
1108}
1109
1110//_______________________________________________________________________________
1111
1112void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1113
1114 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1115 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1116 // has a width of ~5 sigmas. Two band needed for opening angle considerations
76bf81f2 1117 TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1118 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1119 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
931f8b00 1120
1121 Double_t deltaM=candDstar->DeltaInvMass();
1122 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
ed5a88ff 1123 Double_t z=Z(candDstar,jet);
931f8b00 1124 Double_t ptD=candDstar->Pt();
b2705b43 1125
931f8b00 1126 Double_t dPhi=jet->Phi()-candDstar->Phi();
b2705b43 1127
1128 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1129 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
931f8b00 1130
ed5a88ff 1131 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
2942f542 1132 Double_t point[7]={z,dPhi,fPtJet,ptD,deltaM,static_cast<Double_t>(isSideBand), static_cast<Double_t>(fIsDInJet ? 1 : 0)};
931f8b00 1133 //fill the background histogram with the side bands or when looking at MC with the real background
1134 if(isSideBand==1){
1135 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1136 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
ed5a88ff 1137 hsDphiz->Fill(point,1.);
b2705b43 1138 if(fIsDInJet){ // evaluate in the near side
1139 hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
931f8b00 1140 }
1141 }
1142}
1143
1144//_______________________________________________________________________________
1145
b2705b43 1146void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg){
931f8b00 1147
1148 //need mass, deltaR, pt jet, ptD
1149 //two cases: D0 and Dstar
1150
1151 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1152 if(!isselected) return;
1153
1154 Double_t ptD=candbg->Pt();
931f8b00 1155
76bf81f2 1156 TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1157 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
931f8b00 1158
1159
1160 if(fCandidateType==kDstartoKpipi){
1161 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1162 Double_t deltaM=dstarbg->DeltaInvMass();
1163 hInvMassptDbg->Fill(deltaM,ptD);
b2705b43 1164 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
931f8b00 1165 }
1166
1167 if(fCandidateType==kD0toKpi){
1168 Double_t masses[2]={0.,0.};
1169 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1170 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1171
1172 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1173 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1174
1175
1176 if(isselected==1 || isselected==3) {
b2705b43 1177 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
931f8b00 1178 hInvMassptDbg->Fill(masses[0],ptD);
1179 }
1180 if(isselected>=2) {
b2705b43 1181 if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
931f8b00 1182 hInvMassptDbg->Fill(masses[1],ptD);
1183 }
1184
1185
1186 }
1187}
1188
1189//_______________________________________________________________________________
1190
1191Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1192 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1193
1194 if(!p1 || !p2) return -1;
1195 Double_t phi1=p1->Phi(),eta1=p1->Eta();
1196 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1197
1198 Double_t dPhi=phi1-phi2;
1199 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1200 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1201
1202 Double_t dEta=eta1-eta2;
1203 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1204 return deltaR;
1205
1206}
1207
1208//_______________________________________________________________________________
1209
1210Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1211
1212 Double_t ptD=candDstar->Pt();
1213 Int_t bin = fCuts->PtBin(ptD);
1214 if (bin < 0){
1215 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1216 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1217 return -1; // out of bounds
1218 }
1219
1220 Double_t invM=candDstar->InvMassD0();
1221 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1222
1223 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1224 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1225
1226 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1227 else return 0;
1228
1229}
b2705b43 1230
1231//_______________________________________________________________________________
1232
1233Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1234 //returns true/false and the array of particles not found among jet constituents
1235
1236 TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1237 TH1I* hIDddaugh =(TH1I*)fOutput->FindObject("hIDddaugh");
1238 TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1239 TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1240
1241 Int_t daughtersID[3];
1242 Int_t ndaugh=0;
1243 Int_t dmatchedID[3]={0,0,0};
1244 Int_t countmatches=0;
1245 if (fCandidateType==kDstartoKpipi){
1246 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1247 AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1248 daughtersID[0]=dzero->GetProngID(0);
1249 daughtersID[1]=dzero->GetProngID(1);
1250 daughtersID[2]=dstar->GetBachelor()->GetID();
1251 ndaugh=3;
1252
1253 dtrks=new AliAODTrack*[3];
1254 dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1255 dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1256 dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1257
1258 //check
1259 if(fillH){
1260 if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID()) hControlDInJ->Fill(6);
1261
1262 hIDddaugh->Fill(daughtersID[0]);
1263 hIDddaugh->Fill(daughtersID[1]);
1264 hIDddaugh->Fill(daughtersID[2]);
1265
1266 }
1267 //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1268 }
1269
1270 if (fCandidateType==kD0toKpi){
1271 daughtersID[0]=charm->GetProngID(0);
1272 daughtersID[1]=charm->GetProngID(1);
1273 ndaugh=2;
1274 if(fillH){
1275 hIDddaugh->Fill(daughtersID[0]);
1276 hIDddaugh->Fill(daughtersID[1]);
1277 }
1278 dtrks=new AliAODTrack*[2];
1279 dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1280 dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1281
1282 }
1283
1284 const Int_t cndaugh=ndaugh;
1285 daughOutOfJetID=new Int_t[cndaugh];
1286
1287 Int_t nchtrk=thejet->GetNumberOfTracks();
1288 for(Int_t itk=0;itk<nchtrk;itk++){
1289 AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1290 Int_t idtkjet=tkjet->GetID();
1291 if(fillH) hIDjetTracks->Fill(idtkjet);
1292 for(Int_t id=0;id<ndaugh;id++){
1293 if(idtkjet==daughtersID[id]) {
1294 dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1295 countmatches++;
1296
1297 }
1298
1299 if(countmatches==ndaugh) break;
1300 }
1301 }
1302 //reverse: include in the array the ID of daughters not matching
1303
1304 for(Int_t id=0;id<ndaugh;id++){
1305 if(dmatchedID[id]==0) {
1306 daughOutOfJetID[id]=daughtersID[id];
1307 if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1308 }
1309 else daughOutOfJetID[id]=0;
1310 }
1311 if(fillH){
1312 if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1313 if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1314 if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1315 }
1316 if(ndaugh!=countmatches){
1317 return kFALSE;
1318 }
1319
1320 return kTRUE;
1321}
1322
1323//_______________________________________________________________________________
1324
1325Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1326
1327 //check the conditions type:
1328 //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1329 //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1330 //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1331
1332 TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
09fce7b3 1333 TH1F* hDRdaughOut=(TH1F*)fOutput->FindObject("hDRdaughOut");
1334
1335 fPmissing[0]=0;
1336 fPmissing[1]=0;
1337 fPmissing[2]=0;
b2705b43 1338
1339 Bool_t testDeltaR=kFALSE;
1340 if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1341
1342 Int_t* daughOutOfJet;
1343 AliAODTrack** charmDaugh;
1344 Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
09fce7b3 1345 if(testDaugh && testDeltaR) {
1346 //Z of candidates with daughters in the jet
1347 ((TH1F*)fOutput->FindObject("hzDinjet"))->Fill(Z(charm,thejet));
1348
1349 }
b2705b43 1350 if(!testDaugh && testDeltaR && fTypeDInJet==2){
1351
1352 Int_t ndaugh=3;
1353 if(fCandidateType==kD0toKpi) ndaugh=2;
1354 Int_t nOut=ndaugh;
1355
b2705b43 1356 for(Int_t id=0;id<ndaugh;id++){
1357 if(daughOutOfJet[id]!=0){ //non-matched daughter
1358 //get track and its momentum
1359 nOut--;
1360 if(fillH) {
1361 hControlDInJ->Fill(3);
1362 if(id<2) hControlDInJ->Fill(4);
1363 if(id==2)hControlDInJ->Fill(5);
09fce7b3 1364 hDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
b2705b43 1365 }
1366 fPmissing[0]+=charmDaugh[id]->Px();
1367 fPmissing[1]+=charmDaugh[id]->Py();
1368 fPmissing[2]+=charmDaugh[id]->Pz();
1369 }
1370
1371 }
1372
1373 //now the D in within the jet
1374 testDaugh=kTRUE;
1375
1376 }
1377
1378 delete[] charmDaugh;
1379
1380 Bool_t result=0;
1381 switch(fTypeDInJet){
1382 case 0:
1383 result=testDeltaR;
09fce7b3 1384 break;
b2705b43 1385 case 1:
1386 result=testDeltaR && testDaugh;
09fce7b3 1387 break;
b2705b43 1388 case 2:
1389 result=testDeltaR && testDaugh;
09fce7b3 1390 break;
1391 default:
1392 AliInfo("Selection type not specified, use 1");
1393 result=testDeltaR && testDaugh;
1394 break;
b2705b43 1395 }
1396 return result;
1397}