]> 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
Merge branch 'feature-movesplit'
[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>
ad6abcae 33#include <TSystem.h>
34#include <TObjectTable.h>
931f8b00 35
36#include "AliAnalysisTaskFlavourJetCorrelations.h"
37#include "AliAODMCHeader.h"
38#include "AliAODHandler.h"
39#include "AliAnalysisManager.h"
40#include "AliLog.h"
41#include "AliEmcalJet.h"
42#include "AliJetContainer.h"
43#include "AliAODRecoDecay.h"
44#include "AliAODRecoCascadeHF.h"
45#include "AliAODRecoDecayHF2Prong.h"
46#include "AliESDtrack.h"
47#include "AliAODMCParticle.h"
48#include "AliPicoTrack.h"
49#include "AliRDHFCutsD0toKpi.h"
50#include "AliRDHFCutsDStartoKpipi.h"
60e4f65e 51#include "AliRhoParameter.h"
cc5dbb3b 52#include "AliParticleContainer.h"
931f8b00 53
54ClassImp(AliAnalysisTaskFlavourJetCorrelations)
55
56
57//_______________________________________________________________________________
58
59AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
0455d618 60AliAnalysisTaskEmcalJet("",kTRUE),
931f8b00 61fUseMCInfo(kTRUE),
62fUseReco(kTRUE),
63fCandidateType(),
64fPDGmother(),
65fNProngs(),
66fPDGdaughters(),
67fBranchName(),
931f8b00 68fCuts(0),
69fMinMass(),
70fMaxMass(),
71fJetArrName(0),
cc5dbb3b 72fTrackArrName(0),
931f8b00 73fCandArrName(0),
74fLeadingJetOnly(kFALSE),
bbb94467 75fJetRadius(0),
76fCandidateArray(0),
0455d618 77fSideBandArray(0),
b2705b43 78fJetOnlyMode(0),
79fPmissing(),
80fPtJet(0),
81fIsDInJet(0),
82fTypeDInJet(2),
ad6abcae 83fTrackArr(0),
84fSwitchOnSB(0),
85fSwitchOnPhiAxis(0),
86fSwitchOnOutOfConeAxis(0),
87fSwitchOnSparses(1),
0b7f0a99 88fNAxesBigSparse(9),
cc5dbb3b 89fJetCont(0),
90fTrackCont(0),
91fClusterCont(0),
0b7f0a99 92fhstat(),
93fhPtJetTrks(),
94fhPhiJetTrks(),
95fhEtaJetTrks(),
96fhEjetTrks(),
97fhPtJet(),
98fhPhiJet(),
99fhEtaJet(),
100fhEjet(),
101fhNtrArr(),
102fhNJetPerEv(),
103fhdeltaRJetTracks(),
104fhsJet(),
e278882f 105fhImpPar(),
0b7f0a99 106fhNDPerEvNoJet(),
107fhptDPerEvNoJet(),
108fhNJetPerEvNoD(),
109fhPtJetPerEvNoD(),
e278882f 110fhVtxX(),
111fhVtxY(),
112fhVtxZ(),
0b7f0a99 113fhsDstandalone(),
114fhInvMassptD(),
115fhDiffSideBand(),
116fhInvMassptDbg(),
117fhPtPion(),
118fhztracksinjet(),
119fhDiffPtTrPtJzNok(),
120fhDiffPtTrPtJzok(),
121fhDiffPzTrPtJzok(),
122fhDiffPzTrPtJzNok(),
123fhNtrkjzNok(),
124fhztracksinjetT(),
125fhControlDInJ(),
126fhIDddaugh(),
127fhIDddaughOut(),
128fhIDjetTracks(),
129fhDRdaughOut(),
130fhzDinjet(),
131fhmissingp(),
132fhMissPi(),
133fhDeltaPtJet(),
134fhRelDeltaPtJet(),
135fhzDT(),
136fhDeltaRD(),
137fhDeltaRptDptj(),
138fhDeltaRptDptjB(),
139fhsDphiz()
e278882f 140
931f8b00 141{
142 //
143 // Default ctor
144 //
0b7f0a99 145 //SetMakeGeneralHistograms(kTRUE)(),
931f8b00 146
147}
148
149//_______________________________________________________________________________
150
a8b2b672 151AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, Bool_t jetOnly) :
0455d618 152AliAnalysisTaskEmcalJet(name,kTRUE),
931f8b00 153fUseMCInfo(kTRUE),
154fUseReco(kTRUE),
155fCandidateType(),
156fPDGmother(),
157fNProngs(),
158fPDGdaughters(),
159fBranchName(),
931f8b00 160fCuts(0),
161fMinMass(),
162fMaxMass(),
163fJetArrName(0),
cc5dbb3b 164fTrackArrName(0),
931f8b00 165fCandArrName(0),
166fLeadingJetOnly(kFALSE),
bbb94467 167fJetRadius(0),
168fCandidateArray(0),
0455d618 169fSideBandArray(0),
b2705b43 170fJetOnlyMode(jetOnly),
171fPmissing(),
172fPtJet(0),
173fIsDInJet(0),
174fTypeDInJet(2),
ad6abcae 175fTrackArr(0),
176fSwitchOnSB(0),
177fSwitchOnPhiAxis(0),
178fSwitchOnOutOfConeAxis(0),
179fSwitchOnSparses(1),
0b7f0a99 180fNAxesBigSparse(9),
cc5dbb3b 181fJetCont(0),
182fTrackCont(0),
183fClusterCont(0),
0b7f0a99 184fhstat(),
185fhPtJetTrks(),
186fhPhiJetTrks(),
187fhEtaJetTrks(),
188fhEjetTrks(),
189fhPtJet(),
190fhPhiJet(),
191fhEtaJet(),
192fhEjet(),
193fhNtrArr(),
194fhNJetPerEv(),
195fhdeltaRJetTracks(),
196fhsJet(),
e278882f 197fhImpPar(),
0b7f0a99 198fhNDPerEvNoJet(),
199fhptDPerEvNoJet(),
200fhNJetPerEvNoD(),
201fhPtJetPerEvNoD(),
e278882f 202fhVtxX(),
203fhVtxY(),
204fhVtxZ(),
0b7f0a99 205fhsDstandalone(),
206fhInvMassptD(),
207fhDiffSideBand(),
208fhInvMassptDbg(),
209fhPtPion(),
210fhztracksinjet(),
211fhDiffPtTrPtJzNok(),
212fhDiffPtTrPtJzok(),
213fhDiffPzTrPtJzok(),
214fhDiffPzTrPtJzNok(),
215fhNtrkjzNok(),
216fhztracksinjetT(),
217fhControlDInJ(),
218fhIDddaugh(),
219fhIDddaughOut(),
220fhIDjetTracks(),
221fhDRdaughOut(),
222fhzDinjet(),
223fhmissingp(),
224fhMissPi(),
225fhDeltaPtJet(),
226fhRelDeltaPtJet(),
227fhzDT(),
228fhDeltaRD(),
229fhDeltaRptDptj(),
230fhDeltaRptDptjB(),
231fhsDphiz()
931f8b00 232{
233 //
234 // Constructor. Initialization of Inputs and Outputs
235 //
236
237 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
238 fCuts=cuts;
239 fCandidateType=candtype;
240 const Int_t nptbins=fCuts->GetNPtBins();
241 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};
242
243 switch(fCandidateType){
244 case 0:
245 fPDGmother=421;
246 fNProngs=2;
247 fPDGdaughters[0]=211;//pi
248 fPDGdaughters[1]=321;//K
249 fPDGdaughters[2]=0; //empty
250 fPDGdaughters[3]=0; //empty
251 fBranchName="D0toKpi";
252 fCandArrName="D0";
253 break;
254 case 1:
255 fPDGmother=413;
256 fNProngs=3;
257 fPDGdaughters[1]=211;//pi soft
258 fPDGdaughters[0]=421;//D0
259 fPDGdaughters[2]=211;//pi fromD0
260 fPDGdaughters[3]=321; // K from D0
261 fBranchName="Dstar";
262 fCandArrName="Dstar";
263
264 if(nptbins<=13){
265 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
266 } else {
267 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
268 }
269 break;
270 default:
271 printf("%d not accepted!!\n",fCandidateType);
272 break;
273 }
274
275 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
276 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
277
0455d618 278 if(fJetOnlyMode){
279 DefineOutput(1,TList::Class()); //histos with jet info
280 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
281 }
282 else{
283 DefineInput(1, TClonesArray::Class());
284 DefineInput(2, TClonesArray::Class());
285
286 DefineOutput(1,TList::Class()); // histos
287 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
288 }
931f8b00 289}
290
291//_______________________________________________________________________________
292
293AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
294 //
295 // destructor
296 //
297
298 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
299
931f8b00 300 delete fCuts;
301
302}
303
304//_______________________________________________________________________________
305
306void AliAnalysisTaskFlavourJetCorrelations::Init(){
307 //
308 // Initialization
309 //
310
311 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
0455d618 312
931f8b00 313 switch(fCandidateType){
314 case 0:
315 {
316 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
317 copyfCuts->SetName("AnalysisCutsDzero");
318 // Post the data
319 PostData(2,copyfCuts);
320 }
321 break;
322 case 1:
323 {
324 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
325 copyfCuts->SetName("AnalysisCutsDStar");
326 // Post the cuts
327 PostData(2,copyfCuts);
328 }
329 break;
330 default:
331 return;
332 }
333
334 return;
335}
336
337//_______________________________________________________________________________
338
339void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
340 // output
341 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
76bf81f2 342 AliAnalysisTaskEmcal::UserCreateOutputObjects();
cc5dbb3b 343
344 fJetCont = GetJetContainer(0);
345 if(fJetCont){
346 fTrackCont = fJetCont->GetParticleContainer();
347 fClusterCont = fJetCont->GetClusterContainer();
348 }
349
350
931f8b00 351 // define histograms
76bf81f2 352 // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
931f8b00 353 DefineHistoForAnalysis();
76bf81f2 354 PostData(1,fOutput);
931f8b00 355
356 return;
357}
358
359//_______________________________________________________________________________
360
b5d0ba13 361Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
931f8b00 362{
76bf81f2 363 // user exec from AliAnalysisTaskEmcal is used
364
931f8b00 365 // Load the event
366 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
367
368 TClonesArray *arrayDStartoD0pi=0;
931f8b00 369
370 if (!aodEvent && AODEvent() && IsStandardAOD()) {
371
372 // In case there is an AOD handler writing a standard AOD, use the AOD
373 // event in memory rather than the input (ESD) event.
374 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
375
376 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
377 // have to taken from the AOD event hold by the AliAODExtension
378 AliAODHandler* aodHandler = (AliAODHandler*)
379 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
380 if(aodHandler->GetExtensions()) {
381 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
382 AliAODEvent *aodFromExt = ext->GetAOD();
383 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
384 }
76473f1c 385 } else if(aodEvent){
931f8b00 386 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
387 }
388
389 if (!arrayDStartoD0pi) {
390 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
391 // return;
392 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
393
394 TClonesArray* mcArray = 0x0;
da94b30f 395 if (fUseMCInfo) { //not used at the moment,uncomment return if you use
931f8b00 396 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
397 if (!mcArray) {
398 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
da94b30f 399 //return kFALSE;
931f8b00 400 }
401 }
402
403 //retrieve jets
cc5dbb3b 404 //this is a duplication of fTrackCont, but is is used in the loop of line 598 and changing it needs a thorough test
405 fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrName));
931f8b00 406 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
cc5dbb3b 407 //fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
408 //fJetContainer=GetJetContainer(0);
409 //if(!fJetContainer) {
410 // AliError("Jet Container 0 not found");
411 // return kFALSE;
412 //}
413 fJetRadius=fJetCont->GetJetRadius();
931f8b00 414
b2705b43 415 if(!fTrackArr){
931f8b00 416 AliInfo(Form("Could not find the track array\n"));
b5d0ba13 417 return kFALSE;
931f8b00 418 }
419
bbb94467 420
c3c1ab31 421 fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
b5d0ba13 422 if (!fCandidateArray) return kFALSE;
da94b30f 423 if ((fCandidateType==1 && fSwitchOnSB) || fUseMCInfo) {
c3c1ab31 424 fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
b5d0ba13 425 if (!fSideBandArray) return kFALSE;
c3c1ab31 426 }
427 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
931f8b00 428
0b7f0a99 429 fhstat->Fill(0);
931f8b00 430
431 // fix for temporary bug in ESDfilter
432 // the AODs with null vertex pointer didn't pass the PhysSel
b5d0ba13 433 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
931f8b00 434
435 //Event selection
436 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
437 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
b5d0ba13 438 if(!iseventselected) return kFALSE;
931f8b00 439
0b7f0a99 440 fhstat->Fill(1);
e278882f 441
442 const AliVVertex *vertex = aodEvent->GetPrimaryVertex();
443 Double_t vtx[3];
444 vertex->GetXYZ(vtx);
445 fhVtxX->Fill(vtx[0]);
446 fhVtxY->Fill(vtx[1]);
447 fhVtxZ->Fill(vtx[2]);
448 //Printf(">>>>>>>>.>>>Primary vertex %.3f,%.3f,%.3f",vtx[0], vtx[1],vtx[2]);
931f8b00 449 //retrieve charm candidates selected
0455d618 450 Int_t candidates = 0;
cc5dbb3b 451 Int_t njets=fJetCont->GetNJets();
1eb2586e 452 //Printf("N jets in this event %d",njets);
0455d618 453 if(!fJetOnlyMode) {
454 candidates = fCandidateArray->GetEntriesFast();
455
456 //trigger on jets
d9126a4c 457 if(njets == 0) {
0b7f0a99 458 fhstat->Fill(6, candidates);
459 fhNDPerEvNoJet->Fill(candidates);
d9126a4c 460 for(Int_t iD=0;iD<candidates;iD++){
bbb94467 461 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
d9126a4c 462 if(!cand) continue;
0b7f0a99 463 fhptDPerEvNoJet->Fill(cand->Pt());
d9126a4c
CB		 464
465 }
b5d0ba13 466 return kFALSE;
d9126a4c
CB		 467
468 }
1353fa73 469
470 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
471
472 for(Int_t ic = 0; ic < candidates; ic++) {
473
474 // D* candidates
bbb94467 475 AliAODRecoDecayHF* charm=0x0;
476 AliAODRecoCascadeHF* dstar=0x0;
477
478
479 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
1353fa73 480 if(!charm) continue;
bbb94467 481
482 if(fCandidateType==kDstartoKpipi){
483 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
484 if(!dstar) continue;
485 }
486
0b7f0a99 487 fhstat->Fill(2);
1353fa73 488
489 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
490
491 if(fCandidateType==kDstartoKpipi) {
8577f7bc 492 if(fUseReco){
493 Double_t deltamass= dstar->DeltaInvMass();
494 candsparse[3]=deltamass;
495 } else candsparse[3] = 0.145; //for generated
0b7f0a99 496 if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
1353fa73 497 }
498 if(fCandidateType==kD0toKpi){
8577f7bc 499 if(fUseReco){
500 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
501 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
502 //not a further selection,just to get the value of isselected for the D0
503 Double_t masses[2];
504 Int_t pdgdaughtersD0[2]={211,321};//pi,K
505 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
506
507 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
508 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
509 if(isselected==1 || isselected==3) {
510 candsparse[3]=masses[0];
0b7f0a99 511 if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
8577f7bc 512 }
513 if(isselected>=2){
514 candsparse[3]=masses[1];
0b7f0a99 515 if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
8577f7bc 516
517 }
518 } else { //generated
519 Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
520 candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
0b7f0a99 521 if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
1353fa73 522 }
523 }
524 }
0455d618 525 }
60e4f65e 526
527 //Background Subtraction for jets
d8324fea 528 //If there's no background subtraction rhoval=0 and momentum is simply not took into account
60e4f65e 529 AliRhoParameter *rho = 0;
530 Double_t rhoval = 0;
4b0942ca 531 if (!fJetCont->GetRhoName().IsNull()) {
532 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
60e4f65e 533 if(rho) rhoval = rho->GetVal();
534 }
535
0455d618 536
931f8b00 537 // we start with jets
538 Double_t ejet = 0;
539 Double_t phiJet = 0;
540 Double_t etaJet = 0;
541 Double_t ptjet = 0;
542 Double_t leadingJet =0;
0455d618 543 Double_t pointJ[6];
931f8b00 544
cc5dbb3b 545 Int_t ntrarr=fTrackCont->GetNParticles();
0b7f0a99 546 fhNtrArr->Fill(ntrarr);
931f8b00 547
548 for(Int_t i=0;i<ntrarr;i++){
e278882f 549 AliAODTrack *jtrack=static_cast<AliAODTrack*>(fTrackCont->GetParticle(i));
931f8b00 550 //reusing histograms
128389ca 551 if(!jtrack || jtrack->IsMuonTrack()) continue; // added check on IsMuonTrack because in some cases the DCA() gives crazy YAtDCA values that issue floating point exception
0b7f0a99 552 fhPtJetTrks->Fill(jtrack->Pt());
553 fhPhiJetTrks->Fill(jtrack->Phi());
554 fhEtaJetTrks->Fill(jtrack->Eta());
555 fhEjetTrks->Fill(jtrack->E());
e278882f 556 Double_t vDCAglobalxy;
557 //Double_t vDCAglobalz;
558 if(jtrack->IsGlobalConstrained()){
559 // retrieve impact parameter
560 Double_t pos[3];
561
562 jtrack->GetXYZ(pos);
563
564 vDCAglobalxy = pos[0] - vtx[0]; //should be impact parameter in transverse plane
565 //vDCAglobalz = pos[1] - vtx[1]; //should be impact parameter in z direction
128389ca 566 } else {
567 vDCAglobalxy=jtrack->DCA();
568 //vDCAglobalz=jtrack->ZAtDCA();
e278882f 569 }
570 fhImpPar->Fill(vDCAglobalxy,jtrack->Pt());
571 //Printf("Track position %.3f,%.3f,%.3f",pos[0],pos[1],pos[2]);
931f8b00 572 }
573
574
575 //option to use only the leading jet
576 if(fLeadingJetOnly){
577 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
cc5dbb3b 578 AliEmcalJet* jetL = (AliEmcalJet*)fJetCont->GetJet(iJetsL);
d8324fea 579 ptjet = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
931f8b00 580 if(ptjet>leadingJet ) leadingJet = ptjet;
581
582 }
583 }
584
585 Int_t cntjet=0;
586 //loop over jets and consider only the leading jet in the event
09fce7b3 587 for (Int_t iJets = 0; iJets<njets; iJets++) {
588 fPmissing[0]=0;
589 fPmissing[1]=0;
590 fPmissing[2]=0;
591
931f8b00 592 //Printf("Jet N %d",iJets);
cc5dbb3b 593 AliEmcalJet* jet = (AliEmcalJet*)fJetCont->GetJet(iJets);
931f8b00 594 //Printf("Corr task Accept Jet");
595 if(!AcceptJet(jet)) {
0b7f0a99 596 fhstat->Fill(5);
931f8b00 597 continue;
598 }
599
600 //jets variables
601 ejet = jet->E();
602 phiJet = jet->Phi();
603 etaJet = jet->Eta();
d8324fea 604 fPtJet = jet->Pt() - jet->Area()*rhoval; //It takes into account the background subtraction
b2705b43 605 Double_t origPtJet=fPtJet;
09fce7b3 606
0455d618 607 pointJ[0] = phiJet;
608 pointJ[1] = etaJet;
d8324fea 609 pointJ[2] = ptjet - jet->Area()*rhoval; //It takes into account the background subtraction
0455d618 610 pointJ[3] = ejet;
611 pointJ[4] = jet->GetNumberOfConstituents();
612 pointJ[5] = jet->Area();
931f8b00 613
614 // choose the leading jet
615 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
616 //used for normalization
0b7f0a99 617 fhstat->Fill(3);
931f8b00 618 cntjet++;
619 // fill energy, eta and phi of the jet
0b7f0a99 620 fhEjet ->Fill(ejet);
621 fhPhiJet ->Fill(phiJet);
622 fhEtaJet ->Fill(etaJet);
623 fhPtJet ->Fill(fPtJet);
624 if(fJetOnlyMode && fSwitchOnSparses) fhsJet->Fill(pointJ,1);
931f8b00 625 //loop on jet particles
09fce7b3 626 Int_t ntrjet= jet->GetNumberOfTracks();
627 Double_t sumPtTracks=0,sumPzTracks=0;
628 Int_t zg1jtrk=0;
931f8b00 629 for(Int_t itrk=0;itrk<ntrjet;itrk++){
630
b2705b43 631 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
0b7f0a99 632 fhdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
09fce7b3 633 Double_t ztrackj=Z(jetTrk,jet);
0b7f0a99 634 fhztracksinjet->Fill(ztrackj);
09fce7b3 635 sumPtTracks+=jetTrk->Pt();
636 sumPzTracks+=jetTrk->Pz();
637 if(ztrackj>1){
638 zg1jtrk++;
639 }
640
641 Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
0b7f0a99 642 fhztracksinjetT->Fill(ztrackjTr);
09fce7b3 643
931f8b00 644 }//end loop on jet tracks
645
0b7f0a99 646 fhNtrkjzNok->Fill(zg1jtrk);
09fce7b3 647 Double_t diffpt=origPtJet-sumPtTracks;
648 Double_t diffpz=jet->Pz()-sumPzTracks;
649 if(zg1jtrk>0){
0b7f0a99 650 fhDiffPtTrPtJzNok->Fill(diffpt);
651 fhDiffPzTrPtJzNok->Fill(diffpz);
09fce7b3 652
653 }else{
0b7f0a99 654 fhDiffPtTrPtJzok->Fill(diffpt);
655 fhDiffPzTrPtJzok->Fill(diffpz);
09fce7b3 656 }
657
d9126a4c 658 if(candidates==0){
5a51a395 659
0b7f0a99 660 fhPtJetPerEvNoD->Fill(fPtJet);
d9126a4c 661 }
0455d618 662 if(!fJetOnlyMode) {
b2705b43 663 //Printf("N candidates %d ", candidates);
664 for(Int_t ic = 0; ic < candidates; ic++) {
0b7f0a99 665 fhstat->Fill(7);
b2705b43 666 // D* candidates
667 AliVParticle* charm=0x0;
668 charm=(AliVParticle*)fCandidateArray->At(ic);
669 if(!charm) continue;
670 AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
671 fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
672 if (fIsDInJet) FlagFlavour(jet);
0b7f0a99 673 if (jet->TestFlavourTag(AliEmcalJet::kDStar) || jet->TestFlavourTag(AliEmcalJet::kD0)) fhstat->Fill(4);
b2705b43 674
675 //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.
676
677 Double_t pjet[3];
678 jet->PxPyPz(pjet);
d8324fea 679 //It corrects the jet momentum if it was asked for jet background subtraction
60e4f65e 680 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
681 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
682 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 683
684 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 685 RecalculateMomentum(pjet,fPmissing);
d8324fea 686 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
4b0942ca 687 if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
09fce7b3 688
689 //debugging histograms
d8324fea 690 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
0b7f0a99 691 for(Int_t k=0;k<3;k++) fhMissPi[k]->Fill(fPmissing[k]);
b2705b43 692
0b7f0a99 693 fhmissingp->Fill(pmissing);
b2705b43 694 Double_t ptdiff=fPtJet-origPtJet;
0b7f0a99 695 fhDeltaPtJet->Fill(ptdiff);
696 fhRelDeltaPtJet->Fill(ptdiff/origPtJet);
b2705b43 697
698 FillHistogramsRecoJetCorr(charm, jet, aodEvent);
699
700 }//end loop on candidates
931f8b00 701
b2705b43 702 //retrieve side band background candidates for Dstar
703 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
931f8b00 704
b2705b43 705 for(Int_t ib=0;ib<nsbcand;ib++){
706 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
707 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
708 if(!sbcand) continue;
ce11164d 709
b2705b43 710 fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
711 Double_t pjet[3];
712 jet->PxPyPz(pjet);
d8324fea 713 //It corrects the jet momentum if it was asked for jet background subtraction
60e4f65e 714 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
715 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
716 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 717
718 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 719 RecalculateMomentum(pjet,fPmissing);
d8324fea 720 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
4b0942ca 721 if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
b2705b43 722 SideBandBackground(sbcand,jet);
723
724 }
725 if(fUseMCInfo){
da94b30f 726
b2705b43 727 AliAODRecoDecayHF* charmbg = 0x0;
a1a0a01b 728 charmbg=(AliAODRecoDecayHF*)fSideBandArray->At(ib);
b2705b43 729 if(!charmbg) continue;
0b7f0a99 730 fhstat->Fill(8);
b2705b43 731 fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
da94b30f 732 if (fIsDInJet) {
733 FlagFlavour(jet); //this are backgroud HF jets, but flagged as signal at the moment. Can use the bkg flavour flag in the future. This info is not stored now a part in the jet
0b7f0a99 734 fhstat->Fill(9);
da94b30f 735 }
b2705b43 736 Double_t pjet[3];
737 jet->PxPyPz(pjet);
d8324fea 738 //It corrects the jet momentum if it was asked for jet background subtraction
a1a0a01b 739 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
740 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
741 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 742
743 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 744 RecalculateMomentum(pjet,fPmissing);
d8324fea 745 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
4b0942ca 746 if((jet->Pt()-jet->Area()*rhoval)<0) fPtJet = -fPtJet;
a1a0a01b 747 MCBackground(charmbg,jet);
b2705b43 748 }
931f8b00 749 }
750 }
751 } // end of jet loop
752
0b7f0a99 753 fhNJetPerEv->Fill(cntjet);
754 if(candidates==0) fhNJetPerEvNoD->Fill(cntjet);
76bf81f2 755 PostData(1,fOutput);
b5d0ba13 756 return kTRUE;
931f8b00 757}
758
759//_______________________________________________________________________________
760
761void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
762{
763 // The Terminate() function is the last function to be called during
764 // a query. It always runs on the client, it can be used to present
765 // the results graphically or save the results to file.
766
767 Info("Terminate"," terminate");
768 AliAnalysisTaskSE::Terminate();
769
76bf81f2 770 fOutput = dynamic_cast<TList*> (GetOutputData(1));
771 if (!fOutput) {
772 printf("ERROR: fOutput not available\n");
931f8b00 773 return;
774 }
775}
776
777//_______________________________________________________________________________
778
779void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
780 Float_t mass=0;
781 Int_t abspdg=TMath::Abs(pdg);
782
783 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
784 // compute the Delta mass for the D*
785 if(fCandidateType==kDstartoKpipi){
786 Float_t mass1=0;
787 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
788 mass = mass-mass1;
789 }
790
791 fMinMass = mass-range;
792 fMaxMass = mass+range;
793
794 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
795 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
796}
797
798//_______________________________________________________________________________
799
800void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
801 if(uplimit>lowlimit)
802 {
803 fMinMass = lowlimit;
804 fMaxMass = uplimit;
805 }
806 else{
807 printf("Error! Lower limit larger than upper limit!\n");
808 fMinMass = uplimit - uplimit*0.2;
809 fMaxMass = uplimit;
810 }
811}
812
813//_______________________________________________________________________________
814
815Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
816 if(nptbins>30) {
817 AliInfo("Maximum number of bins allowed is 30!");
818 return kFALSE;
819 }
820 if(!width) return kFALSE;
821 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
822 return kTRUE;
823}
824
825//_______________________________________________________________________________
826
09fce7b3 827Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
931f8b00 828 if(!part || !jet){
829 printf("Particle or jet do not exist!\n");
830 return -99;
831 }
832 Double_t p[3],pj[3];
833 Bool_t okpp=part->PxPyPz(p);
834 Bool_t okpj=jet->PxPyPz(pj);
60e4f65e 835
836 //Background Subtraction
d8324fea 837 //It corrects the each component of the jet momentum for Z calculation
60e4f65e 838 AliRhoParameter *rho = 0;
839 Double_t rhoval = 0;
4b0942ca 840 if (!fJetCont->GetRhoName().IsNull()) {
841 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
60e4f65e 842 if(rho){
843 rhoval = rho->GetVal();
844 pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
845 pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
846 pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
847 }
848 }
849
850
851
931f8b00 852 if(!okpp || !okpj){
853 printf("Problems getting momenta\n");
854 return -999;
855 }
b2705b43 856
857 RecalculateMomentum(pj, fPmissing);
09fce7b3 858 if(transverse) return ZT(p,pj);
859 else
860 return Z(p,pj);
861}
862
863//_______________________________________________________________________________
864Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
865
b2705b43 866 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
867 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
931f8b00 868 return z;
869}
870
09fce7b3 871
872//_______________________________________________________________________________
873Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
874
875 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
876 Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
877 return z;
878}
879
931f8b00 880//_______________________________________________________________________________
881
b2705b43 882void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
883
884 pj[0]+=pmissing[0];
885 pj[1]+=pmissing[1];
886 pj[2]+=pmissing[2];
887
888}
889
890//_______________________________________________________________________________
891
931f8b00 892Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
893
894 // Statistics
da94b30f 895 Int_t nbins=8;
896 if(fUseMCInfo) nbins+=2;
897
0b7f0a99 898 fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
899 fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
900 fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
901 fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
902 fhstat->GetXaxis()->SetBinLabel(4,"N jets");
903 fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
904 fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
905 fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
906 fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
da94b30f 907 if(fUseMCInfo) {
0b7f0a99 908 fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
909 fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
910 fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
911 fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
da94b30f 912 }
0b7f0a99 913 fhstat->SetNdivisions(1);
914 fOutput->Add(fhstat);
931f8b00 915
e278882f 916 const Int_t nbinsmass=300;
4b0942ca 917 const Int_t nbinspttrack=500;
918 Int_t nbinsptjet=500;
919 if(!fJetCont->GetRhoName().IsNull()) nbinsptjet=400;
931f8b00 920 const Int_t nbinsptD=100;
921 const Int_t nbinsz=100;
922 const Int_t nbinsphi=200;
1353fa73 923 const Int_t nbinseta=100;
ad6abcae 924
925 //binning for THnSparse
926 const Int_t nbinsSpsmass=50;
927 const Int_t nbinsSpsptjet=100;
928 const Int_t nbinsSpsptD=50;
929 const Int_t nbinsSpsz=100;
930 const Int_t nbinsSpsphi=100;
931 const Int_t nbinsSpseta=60;
932 const Int_t nbinsSpsContrib=100;
933 const Int_t nbinsSpsA=100;
934
4b0942ca 935 const Float_t pttracklims[2]={0.,200.};
936 Float_t ptjetlims[2]={0.,200.};
937 if(!fJetCont->GetRhoName().IsNull()) ptjetlims[0]=-200.;
931f8b00 938 const Float_t ptDlims[2]={0.,50.};
939 const Float_t zlims[2]={0.,1.2};
940 const Float_t philims[2]={0.,6.3};
1353fa73 941 const Float_t etalims[2]={-1.5,1.5};
0455d618 942 const Int_t nContriblims[2]={0,100};
943 const Float_t arealims[2]={0.,2};
1353fa73 944
931f8b00 945 // jet related fistograms
946
0b7f0a99 947 fhEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
948 fhEjetTrks->Sumw2();
949 fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
950 fhPhiJetTrks->Sumw2();
951 fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
952 fhEtaJetTrks->Sumw2();
4b0942ca 953 fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
0b7f0a99 954 fhPtJetTrks->Sumw2();
955
956 fhEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
957 fhEjet->Sumw2();
958 fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
959 fhPhiJet->Sumw2();
960 fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
961 fhEtaJet->Sumw2();
5dc81ccb 962 fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
0b7f0a99 963 fhPtJet->Sumw2();
964
965 fhdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
966 fhdeltaRJetTracks->Sumw2();
967
968 fhNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
969 fhNtrArr->Sumw2();
970
971 fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
972 fhNJetPerEv->Sumw2();
973
e278882f 974 fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used
975 //NB at the moment fillet with jet track imp par!!!
976 fhVtxX = new TH1F("hVtxX", "Vertex X;vtx_x",100,-0.5,0.5);
977 fhVtxY = new TH1F("hVtxY", "Vertex Y;vtx_y",100,-0.5,0.5);
978 fhVtxZ = new TH1F("hVtxZ", "Vertex Z;vtx_z",100,-10,10);
979 /*
980 if(fUseMCInfo){
981 //understand how to do this. At the moment fhImpPar is filled with the impact parameter of jet tracks even if it is written differently in the fhImpPar definition
982 fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used
983 fOutput->Add(fhImpParB);
984
985 }
986 */
0b7f0a99 987
988 fOutput->Add(fhEjetTrks);
989 fOutput->Add(fhPhiJetTrks);
990 fOutput->Add(fhEtaJetTrks);
991 fOutput->Add(fhPtJetTrks);
992 fOutput->Add(fhEjet);
993 fOutput->Add(fhPhiJet);
994 fOutput->Add(fhEtaJet);
995 fOutput->Add(fhPtJet);
996 fOutput->Add(fhdeltaRJetTracks);
997 fOutput->Add(fhNtrArr);
998 fOutput->Add(fhNJetPerEv);
e278882f 999 fOutput->Add(fhImpPar);
1000 fOutput->Add(fhVtxX);
1001 fOutput->Add(fhVtxY);
1002 fOutput->Add(fhVtxZ);
ad6abcae 1003 if(fJetOnlyMode && fSwitchOnSparses){
0455d618 1004 //thnsparse for jets
1005 const Int_t nAxis=6;
5dc81ccb 1006 const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
1007 const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
1008 const Double_t
1009 maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
1010 fhsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
1011 fhsJet->Sumw2();
931f8b00 1012
0b7f0a99 1013 fOutput->Add(fhsJet);
0455d618 1014
931f8b00 1015 }
931f8b00 1016
0455d618 1017 if(!fJetOnlyMode){
09fce7b3 1018
1019 //debugging histograms
0b7f0a99 1020 fhControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
1021 fhControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
1022 fhControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
1023 fhControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
1024 fhControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
1025 fhControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
1026 fhControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
1027 fhControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
1028 fhControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
09fce7b3 1029
0b7f0a99 1030 fhControlDInJ->SetNdivisions(1);
1031 fhControlDInJ->GetXaxis()->SetLabelSize(0.05);
1032 fOutput->Add(fhControlDInJ);
09fce7b3 1033
0b7f0a99 1034 fhmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
1035 fOutput->Add(fhmissingp);
09fce7b3 1036
0b7f0a99 1037 fhMissPi=new TH1F*[3];
09fce7b3 1038 for(Int_t k=0;k<3;k++) {
0b7f0a99 1039 fhMissPi[k]=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
1040 fOutput->Add(fhMissPi[k]);
09fce7b3 1041 }
0b7f0a99 1042 fhDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
09fce7b3 1043
0b7f0a99 1044 fOutput->Add(fhDeltaPtJet);
1045 fhRelDeltaPtJet=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.);
1046 fOutput->Add(fhRelDeltaPtJet);
09fce7b3 1047
0b7f0a99 1048 fhzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
1049 fOutput->Add(fhzDinjet);
09fce7b3 1050 //frag func of tracks in the jet
0b7f0a99 1051 fhztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
1052 fOutput->Add(fhztracksinjet);
09fce7b3 1053
1054 //check jet and tracks
0b7f0a99 1055 fhDiffPtTrPtJzok=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);
1056 fOutput->Add(fhDiffPtTrPtJzok);
1057 fhDiffPtTrPtJzNok=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);
1058 fOutput->Add(fhDiffPtTrPtJzNok);
1059 fhDiffPzTrPtJzok=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);
1060 fOutput->Add(fhDiffPzTrPtJzok);
1061 fhDiffPzTrPtJzNok=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);
1062 fOutput->Add(fhDiffPzTrPtJzNok);
1063 fhNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
1064 fOutput->Add(fhNtrkjzNok);
09fce7b3 1065
1066 //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
0b7f0a99 1067 fhzDT=new TH1F("hzDT", Form("Z of D %s in jet in transverse components;(p_{T}^{D} dot p_{T}^{jet})/p_{T}^{jet}^{2} ", fUseMCInfo ? "(S+B)" : ""),nbinsz,zlims[0],zlims[1]);
1068 fOutput->Add(fhzDT);
1069 fhztracksinjetT=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]);
1070 fOutput->Add(fhztracksinjetT);
09fce7b3 1071
0b7f0a99 1072 fhIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
1073 fOutput->Add(fhIDddaugh);
1074 fhIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
1075 fOutput->Add(fhIDddaughOut);
1076 fhIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
1077 fOutput->Add(fhIDjetTracks);
09fce7b3 1078
0b7f0a99 1079 fhDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
1080 fOutput->Add(fhDRdaughOut);
09fce7b3 1081
1082
0455d618 1083 if(fCandidateType==kDstartoKpipi)
1084 {
da94b30f 1085 if(fSwitchOnSB){
0b7f0a99 1086 fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1087 fhDiffSideBand->SetStats(kTRUE);
1088 fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
1089 fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1090 fhDiffSideBand->Sumw2();
1091 fOutput->Add(fhDiffSideBand);
da94b30f 1092 }
0455d618 1093
0b7f0a99 1094 fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
1095 fhPtPion->SetStats(kTRUE);
1096 fhPtPion->GetXaxis()->SetTitle("GeV/c");
1097 fhPtPion->GetYaxis()->SetTitle("Entries");
1098 fhPtPion->Sumw2();
1099 fOutput->Add(fhPtPion);
0455d618 1100
1101 }
1102 // D related histograms
0b7f0a99 1103 fhNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
1104 fhNDPerEvNoJet->Sumw2();
1105 fOutput->Add(fhNDPerEvNoJet);
0455d618 1106
0b7f0a99 1107 fhptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
1108 fhptDPerEvNoJet->Sumw2();
1109 fOutput->Add(fhptDPerEvNoJet);
0455d618 1110
ad6abcae 1111 if(fSwitchOnSparses){
1112 const Int_t nAxisD=4;
1113 const Int_t nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
1114 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
1115 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
0b7f0a99 1116 fhsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
1117 fhsDstandalone->Sumw2();
ad6abcae 1118
0b7f0a99 1119 fOutput->Add(fhsDstandalone);
ad6abcae 1120 }
0455d618 1121
ad6abcae 1122 /*
0455d618 1123 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]);
1124 hPtJetWithD->Sumw2();
1125 //for the MC this histogram is filled with the real background
1126 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]);
1127 hPtJetWithDsb->Sumw2();
1128
ad6abcae 1129 fOutput->Add(hPtJetWithD);
1130 fOutput->Add(hPtJetWithDsb);
1131
1132 */
0b7f0a99 1133 fhNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
1134 fhNJetPerEvNoD->Sumw2();
0455d618 1135
0b7f0a99 1136 fhPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1137 fhPtJetPerEvNoD->Sumw2();
0455d618 1138
0b7f0a99 1139 fOutput->Add(fhNJetPerEvNoD);
1140 fOutput->Add(fhPtJetPerEvNoD);
0455d618 1141
0b7f0a99 1142 fhDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1143 fhDeltaRD->Sumw2();
1144 fOutput->Add(fhDeltaRD);
0455d618 1145
0b7f0a99 1146 fhDeltaRptDptj=new TH3F("hDeltaRptDptj",Form("#Delta R distribution of D candidates %s selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)", fUseMCInfo ? "(S)" : ""),100, 0.,5.,nbinsptD, ptDlims[0],ptDlims[1],nbinsptjet,ptjetlims[0],ptjetlims[1]);
1147 fhDeltaRptDptj->Sumw2();
1148 fOutput->Add(fhDeltaRptDptj);
5a51a395 1149
1150 if(fUseMCInfo){
0b7f0a99 1151 fhDeltaRptDptjB=new TH3F("hDeltaRptDptjB",Form("#Delta R distribution of D candidates (B) selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)"),100, 0.,5.,nbinsptD, ptDlims[0],ptDlims[1],nbinsptjet,ptjetlims[0],ptjetlims[1]);
1152 fhDeltaRptDptjB->Sumw2();
1153 fOutput->Add(fhDeltaRptDptjB);
5a51a395 1154 }
1155
0455d618 1156 //background (side bands for the Dstar and like sign for D0)
cc5dbb3b 1157 AliJetContainer *jetCont=GetJetContainer(0);
1158 if(!jetCont){
1159 Printf("Container 0 not found, try with name %s", fJetArrName.Data());
1160 jetCont=GetJetContainer(fJetArrName);
1161 if(!jetCont) Printf("NOT FOUND AGAIN");
1162 return kFALSE;
1163 }
1164 Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
1165 //fJetRadius=jetCont->GetJetRadius();
0b7f0a99 1166 fhInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
1167 fhInvMassptD->SetStats(kTRUE);
1168 fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1169 fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1170 fhInvMassptD->Sumw2();
0455d618 1171
0b7f0a99 1172 fOutput->Add(fhInvMassptD);
0455d618 1173
1174 if(fUseMCInfo){
0b7f0a99 1175 fhInvMassptDbg = 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]);
1176 fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1177 fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1178 fhInvMassptDbg->Sumw2();
1179 fOutput->Add(fhInvMassptDbg);
0455d618 1180
1181 }
0455d618 1182
ad6abcae 1183 if(fSwitchOnSparses){
1184 Double_t pi=TMath::Pi(), philims2[2];
1185 philims2[0]=-pi*1./2.;
1186 philims2[1]=pi*3./2.;
0b7f0a99 1187 fhsDphiz=0x0; //definition below according to the switches
ad6abcae 1188
1189 if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1190 AliInfo("Creating a 9 axes container: This might cause large memory usage");
1191 const Int_t nAxis=9;
5dc81ccb 1192 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1193 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1194 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
ad6abcae 1195 fNAxesBigSparse=nAxis;
1196
0b7f0a99 1197 fhsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. SB? D within the jet cone?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
ad6abcae 1198 }
1199
1200 if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1201 fSwitchOnPhiAxis=0;
1202 fSwitchOnOutOfConeAxis=0;
1203 fSwitchOnSB=0;
a1a0a01b 1204 if(fUseMCInfo){
1205 AliInfo("Creating a 7 axes container (MB background candidates)");
da94b30f 1206 const Int_t nAxis=7;
5dc81ccb 1207 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1208 const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1209 const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
a1a0a01b 1210 fNAxesBigSparse=nAxis;
0b7f0a99 1211 fhsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
a1a0a01b 1212
1213 } else{
1214 AliInfo("Creating a 6 axes container");
1215 const Int_t nAxis=6;
5dc81ccb 1216 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1217 const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1218 const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
a1a0a01b 1219 fNAxesBigSparse=nAxis;
1220
0b7f0a99 1221 fhsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass., D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
a1a0a01b 1222 }
ad6abcae 1223 }
0b7f0a99 1224 if(!fhsDphiz) AliFatal("No THnSparse created");
1225 fhsDphiz->Sumw2();
ad6abcae 1226
0b7f0a99 1227 fOutput->Add(fhsDphiz);
ad6abcae 1228 }
0455d618 1229 }
76bf81f2 1230 PostData(1, fOutput);
931f8b00 1231
1232 return kTRUE;
1233}
1234
1235//_______________________________________________________________________________
1236
8577f7bc 1237void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet, AliAODEvent* aodEvent){
931f8b00 1238
1239 Double_t ptD=candidate->Pt();
d8324fea 1240 Double_t deltaR=DeltaR(jet,candidate);
b2705b43 1241 Double_t phiD=candidate->Phi();
1242 Double_t deltaphi = jet->Phi()-phiD;
1243 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1244 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
931f8b00 1245 Double_t z=Z(candidate,jet);
09fce7b3 1246 /*
1247 if(z>1) {
1248 ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1249 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1250
1251 for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1252
1253 ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1254 Double_t ptdiff=fPtJet-jet->Pt();
1255 ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1256 ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1257
1258
1259 }
1260 */
0b7f0a99 1261 if(fIsDInJet)fhzDT->Fill(Z(candidate,jet,kTRUE));
1262
1263 fhDeltaRD->Fill(deltaR);
1264 fhDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
5a51a395 1265
ce11164d 1266 Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1267 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
931f8b00 1268 if(fUseReco){
1269 if(fCandidateType==kD0toKpi) {
1270 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
8577f7bc 1271
ad6abcae 1272 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
931f8b00 1273
1274 }
1275
1276 if(fCandidateType==kDstartoKpipi) {
1277 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
ad6abcae 1278 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
931f8b00 1279
1280 }
1281 } else{ //generated level
1282 //AliInfo("Non reco");
ad6abcae 1283 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
931f8b00 1284 }
1285}
1286
1287//_______________________________________________________________________________
1288
ad6abcae 1289void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent){
931f8b00 1290
931f8b00 1291
1292 Double_t masses[2]={0.,0.};
1293 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1294 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1295
1296 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1297 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1298
ad6abcae 1299 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
ad6abcae 1300 Double_t *point=0x0;
1301 if(fNAxesBigSparse==9){
1302 point=new Double_t[9];
1303 point[0]=z;
1304 point[1]=dPhi;
1305 point[2]=ptj;
1306 point[3]=ptD;
1307 point[4]=masses[0];
1308 point[5]=0;
1309 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1310 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1311 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1312 }
1313 if(fNAxesBigSparse==6){
1314 point=new Double_t[6];
1315 point[0]=z;
1316 point[1]=ptj;
1317 point[2]=ptD;
1318 point[3]=masses[0];
1319 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1320 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
a1a0a01b 1321}
1322 if(fNAxesBigSparse==7){
1323 point=new Double_t[7];
1324 point[0]=z;
1325 point[1]=ptj;
1326 point[2]=ptD;
1327 point[3]=masses[0];
1328 point[4]=0;
1329 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1330 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1331
ad6abcae 1332 }
1070b230 1333 if(!point){
1334 AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1335 return;
1336 }
ad6abcae 1337
e278882f 1338 //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
931f8b00 1339 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1340 if(isselected==1 || isselected==3) {
1341
ad6abcae 1342 //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
931f8b00 1343
b2705b43 1344 FillMassHistograms(masses[0], ptD);
0b7f0a99 1345 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1346 }
1347 if(isselected>=2) {
ad6abcae 1348 //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
931f8b00 1349
b2705b43 1350 FillMassHistograms(masses[1], ptD);
ad6abcae 1351 if(fNAxesBigSparse==9) point[4]=masses[1];
a1a0a01b 1352 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
0b7f0a99 1353 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1354 }
ab1bd143 1355 delete[] point;
931f8b00 1356}
1357
1358//_______________________________________________________________________________
1359
ad6abcae 1360void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
931f8b00 1361 //dPhi and z not used at the moment,but will be (re)added
1362
1363 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1364 Double_t deltamass= dstar->DeltaInvMass();
1365 //Double_t massD0= dstar->InvMassD0();
1366
1367
0b7f0a99 1368 fhPtPion->Fill(softpi->Pt());
931f8b00 1369
ad6abcae 1370 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
ed5a88ff 1371 Int_t isSB=0;//IsDzeroSideBand(dstar);
931f8b00 1372
ad6abcae 1373 //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1374 Double_t *point=0x0;
1375 if(fNAxesBigSparse==9){
1376 point=new Double_t[9];
1377 point[0]=z;
1378 point[1]=dPhi;
1379 point[2]=ptj;
1380 point[3]=ptD;
1381 point[4]=deltamass;
1382 point[5]=static_cast<Double_t>(isSB);
1383 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1384 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1385 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1386 }
1387 if(fNAxesBigSparse==6){
1388 point=new Double_t[6];
1389 point[0]=z;
1390 point[1]=ptj;
1391 point[2]=ptD;
1392 point[3]=deltamass;
1393 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1394 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1395 }
a1a0a01b 1396 if(fNAxesBigSparse==7){
1397 point=new Double_t[7];
1398 point[0]=z;
1399 point[1]=ptj;
1400 point[2]=ptD;
1401 point[3]=deltamass;
1402 point[4]=0;
1403 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1404 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1405 }
931f8b00 1406
1070b230 1407 if(!point){
1408 AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1409 return;
1410 }
1411
ad6abcae 1412 //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
931f8b00 1413
b2705b43 1414 FillMassHistograms(deltamass, ptD);
0b7f0a99 1415 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
ab1bd143 1416 delete[] point;
931f8b00 1417}
1418
1419//_______________________________________________________________________________
1420
ad6abcae 1421void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
931f8b00 1422
1423 Double_t pdgmass=0;
a1a0a01b 1424 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1425 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
ad6abcae 1426 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
0b7f0a99 1427
ad6abcae 1428 //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
5a51a395 1429
ad6abcae 1430 Double_t *point=0x0;
1431 if(fNAxesBigSparse==9){
1432 point=new Double_t[9];
1433 point[0]=z;
1434 point[1]=dPhi;
1435 point[2]=ptjet;
1436 point[3]=ptD;
1437 point[4]=pdgmass;
1438 point[5]=0;
1439 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1440 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1441 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1442 }
1443 if(fNAxesBigSparse==6){
1444 point=new Double_t[6];
1445 point[0]=z;
1446 point[1]=ptjet;
1447 point[2]=ptD;
1448 point[3]=pdgmass;
1449 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1450 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1451 }
a1a0a01b 1452 if(fNAxesBigSparse==7){
da94b30f 1453 point=new Double_t[7];
a1a0a01b 1454 point[0]=z;
1455 point[1]=ptjet;
1456 point[2]=ptD;
1457 point[3]=pdgmass;
1458 point[4]=1;
1459 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1460 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1461 }
931f8b00 1462
1070b230 1463 if(!point){
1464 AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1465 return;
1466 }
a1a0a01b 1467
1468
ad6abcae 1469 if(fNAxesBigSparse==9) point[4]=pdgmass;
a1a0a01b 1470 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
0b7f0a99 1471 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
ad6abcae 1472 //if(fIsDInJet) {
1473 // hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1474 //}
ab1bd143 1475 delete[] point;
931f8b00 1476}
1477
1478//_______________________________________________________________________________
1479
b2705b43 1480void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
931f8b00 1481
b2705b43 1482 if(fIsDInJet) {
0b7f0a99 1483 fhInvMassptD->Fill(mass,ptD);
931f8b00 1484 }
1485}
1486
1353fa73 1487//________________________________________________________________________________
1488
b2705b43 1489void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1490
931f8b00 1491 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1492 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
b2705b43 1493 if (fIsDInJet) jet->AddFlavourTag(tag);
931f8b00 1494
1495 return;
1496
1497}
1498
1499//_______________________________________________________________________________
1500
1501void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1502
1503 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1504 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1505 // has a width of ~5 sigmas. Two band needed for opening angle considerations
0b7f0a99 1506
ad6abcae 1507 //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
931f8b00 1508
ce11164d 1509 Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1510 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1511
931f8b00 1512 Double_t deltaM=candDstar->DeltaInvMass();
1513 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
ed5a88ff 1514 Double_t z=Z(candDstar,jet);
931f8b00 1515 Double_t ptD=candDstar->Pt();
b2705b43 1516
931f8b00 1517 Double_t dPhi=jet->Phi()-candDstar->Phi();
b2705b43 1518
1519 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1520 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
931f8b00 1521
ed5a88ff 1522 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
ad6abcae 1523 //if no SB analysis we should not enter here, so no need of checking the number of axes
1524 Double_t point[9]={z,dPhi,fPtJet,ptD,deltaM,static_cast<Double_t>(isSideBand), static_cast<Double_t>(fIsDInJet ? 1 : 0),static_cast<Double_t>(bDInEMCalAcc? 1 : 0),static_cast<Double_t>(bJetInEMCalAcc? 1 : 0)};
931f8b00 1525 //fill the background histogram with the side bands or when looking at MC with the real background
1526 if(isSideBand==1){
0b7f0a99 1527 fhDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
931f8b00 1528 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
0b7f0a99 1529 if(fSwitchOnSparses) fhsDphiz->Fill(point,1.);
ad6abcae 1530 //if(fIsDInJet){ // evaluate in the near side
1531 // hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1532 //}
931f8b00 1533 }
1534}
1535
1536//_______________________________________________________________________________
1537
a1a0a01b 1538void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
931f8b00 1539
1540 //need mass, deltaR, pt jet, ptD
1541 //two cases: D0 and Dstar
1542
1543 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1544 if(!isselected) return;
1545
1546 Double_t ptD=candbg->Pt();
d8324fea 1547 Double_t deltaR=DeltaR(jet,candbg);
a1a0a01b 1548 Double_t phiD=candbg->Phi();
1549 Double_t deltaphi = jet->Phi()-phiD;
1550 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1551 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1552 Double_t z=Z(candbg,jet);
1553
0b7f0a99 1554 if(fIsDInJet) fhzDT->Fill(Z(candbg,jet,kTRUE));
da94b30f 1555
5a51a395 1556
0b7f0a99 1557
1558
1559 fhDeltaRD->Fill(deltaR);
1560 fhDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
da94b30f 1561
a1a0a01b 1562 Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1563 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1564
ad6abcae 1565 //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
ce11164d 1566
a1a0a01b 1567 Double_t *point=0x0;
1568 if(fNAxesBigSparse==9){
1569 point=new Double_t[9];
1570 point[0]=z;
1571 point[1]=deltaphi;
1572 point[2]=fPtJet;
1573 point[3]=ptD;
1574 point[4]=-999; //set below
1575 point[5]=1;
1576 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1577 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1578 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1579 }
1580
1581 if(fNAxesBigSparse==7){
1582 point=new Double_t[7];
1583 point[0]=z;
1584 point[1]=fPtJet;
1585 point[2]=ptD;
1586 point[3]=-999; //set below
1587 point[4]=1;
1588 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1589 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1590 }
1070b230 1591 if(!point){
1592 AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
1593 return;
1594 }
ce11164d 1595
931f8b00 1596 if(fCandidateType==kDstartoKpipi){
1597 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1598 Double_t deltaM=dstarbg->DeltaInvMass();
0b7f0a99 1599 fhInvMassptDbg->Fill(deltaM,ptD);
ad6abcae 1600 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
a1a0a01b 1601 if(fNAxesBigSparse==9) point[4]=deltaM;
1602 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
0b7f0a99 1603 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1604 }
1605
1606 if(fCandidateType==kD0toKpi){
1607 Double_t masses[2]={0.,0.};
1608 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1609 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1610
1611 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1612 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1613
1614
1615 if(isselected==1 || isselected==3) {
ad6abcae 1616 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
0b7f0a99 1617 fhInvMassptDbg->Fill(masses[0],ptD);
a1a0a01b 1618 if(fNAxesBigSparse==9) point[4]=masses[0];
1619 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
0b7f0a99 1620 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
a1a0a01b 1621 }
931f8b00 1622 if(isselected>=2) {
ad6abcae 1623 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
0b7f0a99 1624 fhInvMassptDbg->Fill(masses[1],ptD);
a1a0a01b 1625 if(fNAxesBigSparse==9) point[4]=masses[1];
1626 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
0b7f0a99 1627 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
a1a0a01b 1628
931f8b00 1629 }
1630
1631
1632 }
ab1bd143 1633 delete[] point;
931f8b00 1634}
1635
1636//_______________________________________________________________________________
1637
d8324fea 1638Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
931f8b00 1639 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
d8324fea 1640 //It recalculates the eta-phi values if it was asked for background subtraction of the jets
931f8b00 1641 if(!p1 || !p2) return -1;
d8324fea 1642
1643 Double_t phi1=p1->Phi(),eta1=p1->Eta();
1644
1645 //It subtracts the backgroud of jets if it was asked for it.
4b0942ca 1646
1647 if (!fJetCont->GetRhoName().IsNull()) {
d8324fea 1648 AliRhoParameter *rho = 0;
4b0942ca 1649 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
d8324fea 1650 if(rho){
1651 Double_t pj[3];
1652 Bool_t okpj=p1->PxPyPz(pj);
1653 if(!okpj){
1654 printf("Problems getting momenta\n");
1655 return -999;
1656 }
1657 Double_t rhoval = rho->GetVal();
1658 pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1659 pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1660 pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1661 //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1662 //[0,pi] if py > 0 and
1663 //[pi,2pi] if py < 0
1664 phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1665 if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1666 eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1667 }
1668 }
1669
931f8b00 1670 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1671
1672 Double_t dPhi=phi1-phi2;
1673 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1674 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1675
1676 Double_t dEta=eta1-eta2;
1677 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1678 return deltaR;
1679
1680}
1681
1682//_______________________________________________________________________________
1683
1684Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1685
1686 Double_t ptD=candDstar->Pt();
1687 Int_t bin = fCuts->PtBin(ptD);
1688 if (bin < 0){
1689 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1690 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1691 return -1; // out of bounds
1692 }
1693
1694 Double_t invM=candDstar->InvMassD0();
1695 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1696
1697 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1698 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1699
1700 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1701 else return 0;
1702
1703}
b2705b43 1704
1705//_______________________________________________________________________________
1706
1707Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1708 //returns true/false and the array of particles not found among jet constituents
1709
b2705b43 1710
1711 Int_t daughtersID[3];
1712 Int_t ndaugh=0;
1713 Int_t dmatchedID[3]={0,0,0};
1714 Int_t countmatches=0;
1715 if (fCandidateType==kDstartoKpipi){
1716 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1717 AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1718 daughtersID[0]=dzero->GetProngID(0);
1719 daughtersID[1]=dzero->GetProngID(1);
1720 daughtersID[2]=dstar->GetBachelor()->GetID();
1721 ndaugh=3;
1722
1723 dtrks=new AliAODTrack*[3];
1724 dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1725 dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1726 dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1727
1728 //check
1729 if(fillH){
0b7f0a99 1730 if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID()) fhControlDInJ->Fill(6);
b2705b43 1731
0b7f0a99 1732 fhIDddaugh->Fill(daughtersID[0]);
1733 fhIDddaugh->Fill(daughtersID[1]);
1734 fhIDddaugh->Fill(daughtersID[2]);
b2705b43 1735
1736 }
1737 //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1738 }
1739
1740 if (fCandidateType==kD0toKpi){
1741 daughtersID[0]=charm->GetProngID(0);
1742 daughtersID[1]=charm->GetProngID(1);
1743 ndaugh=2;
1744 if(fillH){
0b7f0a99 1745 fhIDddaugh->Fill(daughtersID[0]);
1746 fhIDddaugh->Fill(daughtersID[1]);
b2705b43 1747 }
1748 dtrks=new AliAODTrack*[2];
1749 dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1750 dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1751
1752 }
1753
1754 const Int_t cndaugh=ndaugh;
1755 daughOutOfJetID=new Int_t[cndaugh];
1756
1757 Int_t nchtrk=thejet->GetNumberOfTracks();
1758 for(Int_t itk=0;itk<nchtrk;itk++){
1759 AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
e278882f 1760 if(!tkjet) continue;
b2705b43 1761 Int_t idtkjet=tkjet->GetID();
0b7f0a99 1762 if(fillH) fhIDjetTracks->Fill(idtkjet);
b2705b43 1763 for(Int_t id=0;id<ndaugh;id++){
1764 if(idtkjet==daughtersID[id]) {
1765 dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1766 countmatches++;
1767
1768 }
1769
1770 if(countmatches==ndaugh) break;
1771 }
1772 }
1773 //reverse: include in the array the ID of daughters not matching
1774
1775 for(Int_t id=0;id<ndaugh;id++){
1776 if(dmatchedID[id]==0) {
1777 daughOutOfJetID[id]=daughtersID[id];
0b7f0a99 1778 if(fillH) fhIDddaughOut->Fill(daughOutOfJetID[id]);
b2705b43 1779 }
1780 else daughOutOfJetID[id]=0;
1781 }
1782 if(fillH){
0b7f0a99 1783 if((ndaugh-countmatches) == 1) fhControlDInJ->Fill(0);
1784 if((ndaugh-countmatches) == 2) fhControlDInJ->Fill(1);
1785 if((ndaugh-countmatches) == 3) fhControlDInJ->Fill(2);
b2705b43 1786 }
1787 if(ndaugh!=countmatches){
1788 return kFALSE;
1789 }
1790
1791 return kTRUE;
1792}
1793
1794//_______________________________________________________________________________
1795
1796Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1797
1798 //check the conditions type:
1799 //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1800 //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1801 //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
e712e4fb 1802 //type 3 (under development) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet usign the pt-scheme
0b7f0a99 1803
09fce7b3 1804 fPmissing[0]=0;
1805 fPmissing[1]=0;
1806 fPmissing[2]=0;
b2705b43 1807
1808 Bool_t testDeltaR=kFALSE;
1809 if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
128389ca 1810 //for type 3 this check should be redone after the modification of the jet axis
b2705b43 1811
1812 Int_t* daughOutOfJet;
1813 AliAODTrack** charmDaugh;
1814 Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
09fce7b3 1815 if(testDaugh && testDeltaR) {
1816 //Z of candidates with daughters in the jet
0b7f0a99 1817 fhzDinjet->Fill(Z(charm,thejet));
09fce7b3 1818
1819 }
128389ca 1820
1821 TVector3 thejetv; //(x=0,y=0,z=0)
1822 thejetv.SetPtEtaPhi(thejet->Pt(),thejet->Eta(),thejet->Phi());
1823 TVector3 newjet(thejetv);
1824 TVector3 correction; //(x=0,y=0,z=0)
1825 TVector3 charmcand; //(x=0,y=0,z=0)
1826 charmcand.SetPtEtaPhi(charm->Pt(),charm->Eta(),charm->Phi());
1827
e712e4fb 1828 if(!testDaugh && testDeltaR && fTypeDInJet>=2){
b2705b43 1829
1830 Int_t ndaugh=3;
1831 if(fCandidateType==kD0toKpi) ndaugh=2;
1832 Int_t nOut=ndaugh;
1833
b2705b43 1834 for(Int_t id=0;id<ndaugh;id++){
1835 if(daughOutOfJet[id]!=0){ //non-matched daughter
1836 //get track and its momentum
1837 nOut--;
1838 if(fillH) {
0b7f0a99 1839 fhControlDInJ->Fill(3);
1840 if(id<2) fhControlDInJ->Fill(4);
1841 if(id==2)fhControlDInJ->Fill(5);
1842 fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
b2705b43 1843 }
e712e4fb 1844 if(fTypeDInJet==2){
128389ca 1845
1846 newjet.SetX(newjet(0)+charmDaugh[id]->Px());
1847 newjet.SetY(newjet(1)+charmDaugh[id]->Py());
1848 newjet.SetZ(newjet(2)+charmDaugh[id]->Pz());
e712e4fb 1849 }
1850 if(fTypeDInJet==3){
128389ca 1851
e712e4fb 1852 Double_t ptdaug = charmDaugh[id]->Pt();
128389ca 1853 Double_t ptjet = newjet.Perp();
e712e4fb 1854 Double_t ptn = ptjet+ptdaug;
1855 Double_t phidaug = charmDaugh[id]->Phi();
128389ca 1856 Double_t phijet = newjet.Phi();
1857 Double_t phin = (phijet*ptjet+phidaug*ptdaug)/(ptjet+ptdaug);
1858
e712e4fb 1859 Double_t etadaug = charmDaugh[id]->Eta();
128389ca 1860 Double_t etajet = newjet.Eta();
1861 Double_t etan = (etajet*ptjet+etadaug*ptdaug)/(ptjet+ptdaug);
e712e4fb 1862
128389ca 1863 newjet.SetPtEtaPhi(ptn,etan,phin);
1864
e712e4fb 1865 }
128389ca 1866 }
b2705b43 1867 }
ab1bd143 1868
128389ca 1869 correction=newjet-thejetv;
1870 fPmissing[0]=correction(0);
1871 fPmissing[1]=correction(1);
1872 fPmissing[2]=correction(2);
1873
1874 //now the D is within the jet
b2705b43 1875 testDaugh=kTRUE;
128389ca 1876 if(fTypeDInJet==3){ //recalc R(D,jet)
1877 Double_t dRnew=newjet.DeltaR(charmcand);
1878 if(dRnew<fJetRadius) testDeltaR=kTRUE;
1879 else testDeltaR=kFALSE;
1880 }
b2705b43 1881 }
ab1bd143 1882 delete[] daughOutOfJet;
b2705b43 1883 delete[] charmDaugh;
1884
1885 Bool_t result=0;
1886 switch(fTypeDInJet){
1887 case 0:
1888 result=testDeltaR;
09fce7b3 1889 break;
b2705b43 1890 case 1:
1891 result=testDeltaR && testDaugh;
09fce7b3 1892 break;
e712e4fb 1893 case 2: //this case defines fPmissing
1894 result=testDeltaR && testDaugh;
09fce7b3 1895 break;
128389ca 1896 case 3: //this case defines fPmissing and recalculate R(jet,D) with the updated jet axis
1897 result=testDeltaR && testDaugh;
1898 break;
1899
09fce7b3 1900 default:
1901 AliInfo("Selection type not specified, use 1");
1902 result=testDeltaR && testDaugh;
1903 break;
b2705b43 1904 }
128389ca 1905 return result;
b2705b43 1906}
ce11164d 1907
1908Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1909 //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1910
1911 Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1912 Bool_t binEMCal=kTRUE;
1913 Double_t phi=vpart->Phi(), eta=vpart->Eta();
1914 if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1915 if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1916 return binEMCal;
1917
1918
1919}