]> 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
method 3 to add pt missing
[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 }
385 } else {
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;
531 TString sname("Rho");
532 if (!sname.IsNull()) {
533 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
534 if(rho) rhoval = rho->GetVal();
535 }
536
0455d618 537
931f8b00 538 // we start with jets
539 Double_t ejet = 0;
540 Double_t phiJet = 0;
541 Double_t etaJet = 0;
542 Double_t ptjet = 0;
543 Double_t leadingJet =0;
0455d618 544 Double_t pointJ[6];
931f8b00 545
cc5dbb3b 546 Int_t ntrarr=fTrackCont->GetNParticles();
0b7f0a99 547 fhNtrArr->Fill(ntrarr);
931f8b00 548
549 for(Int_t i=0;i<ntrarr;i++){
e278882f 550 AliAODTrack *jtrack=static_cast<AliAODTrack*>(fTrackCont->GetParticle(i));
931f8b00 551 //reusing histograms
cc5dbb3b 552 if(!jtrack) continue;
0b7f0a99 553 fhPtJetTrks->Fill(jtrack->Pt());
554 fhPhiJetTrks->Fill(jtrack->Phi());
555 fhEtaJetTrks->Fill(jtrack->Eta());
556 fhEjetTrks->Fill(jtrack->E());
e278882f 557 Double_t vDCAglobalxy;
558 //Double_t vDCAglobalz;
559 if(jtrack->IsGlobalConstrained()){
560 // retrieve impact parameter
561 Double_t pos[3];
562
563 jtrack->GetXYZ(pos);
564
565 vDCAglobalxy = pos[0] - vtx[0]; //should be impact parameter in transverse plane
566 //vDCAglobalz = pos[1] - vtx[1]; //should be impact parameter in z direction
567 }
568 else{
569 vDCAglobalxy=jtrack->DCA();
570 //vDCAglobalz=jtrack->ZAtDCA();
571 }
572 fhImpPar->Fill(vDCAglobalxy,jtrack->Pt());
573 //Printf("Track position %.3f,%.3f,%.3f",pos[0],pos[1],pos[2]);
931f8b00 574 }
575
576
577 //option to use only the leading jet
578 if(fLeadingJetOnly){
579 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
cc5dbb3b 580 AliEmcalJet* jetL = (AliEmcalJet*)fJetCont->GetJet(iJetsL);
d8324fea 581 ptjet = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
931f8b00 582 if(ptjet>leadingJet ) leadingJet = ptjet;
583
584 }
585 }
586
587 Int_t cntjet=0;
588 //loop over jets and consider only the leading jet in the event
09fce7b3 589 for (Int_t iJets = 0; iJets<njets; iJets++) {
590 fPmissing[0]=0;
591 fPmissing[1]=0;
592 fPmissing[2]=0;
593
931f8b00 594 //Printf("Jet N %d",iJets);
cc5dbb3b 595 AliEmcalJet* jet = (AliEmcalJet*)fJetCont->GetJet(iJets);
931f8b00 596 //Printf("Corr task Accept Jet");
597 if(!AcceptJet(jet)) {
0b7f0a99 598 fhstat->Fill(5);
931f8b00 599 continue;
600 }
601
602 //jets variables
603 ejet = jet->E();
604 phiJet = jet->Phi();
605 etaJet = jet->Eta();
d8324fea 606 fPtJet = jet->Pt() - jet->Area()*rhoval; //It takes into account the background subtraction
b2705b43 607 Double_t origPtJet=fPtJet;
09fce7b3 608
0455d618 609 pointJ[0] = phiJet;
610 pointJ[1] = etaJet;
d8324fea 611 pointJ[2] = ptjet - jet->Area()*rhoval; //It takes into account the background subtraction
0455d618 612 pointJ[3] = ejet;
613 pointJ[4] = jet->GetNumberOfConstituents();
614 pointJ[5] = jet->Area();
931f8b00 615
616 // choose the leading jet
617 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
618 //used for normalization
0b7f0a99 619 fhstat->Fill(3);
931f8b00 620 cntjet++;
621 // fill energy, eta and phi of the jet
0b7f0a99 622 fhEjet ->Fill(ejet);
623 fhPhiJet ->Fill(phiJet);
624 fhEtaJet ->Fill(etaJet);
625 fhPtJet ->Fill(fPtJet);
626 if(fJetOnlyMode && fSwitchOnSparses) fhsJet->Fill(pointJ,1);
931f8b00 627 //loop on jet particles
09fce7b3 628 Int_t ntrjet= jet->GetNumberOfTracks();
629 Double_t sumPtTracks=0,sumPzTracks=0;
630 Int_t zg1jtrk=0;
931f8b00 631 for(Int_t itrk=0;itrk<ntrjet;itrk++){
632
b2705b43 633 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
0b7f0a99 634 fhdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
09fce7b3 635 Double_t ztrackj=Z(jetTrk,jet);
0b7f0a99 636 fhztracksinjet->Fill(ztrackj);
09fce7b3 637 sumPtTracks+=jetTrk->Pt();
638 sumPzTracks+=jetTrk->Pz();
639 if(ztrackj>1){
640 zg1jtrk++;
641 }
642
643 Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
0b7f0a99 644 fhztracksinjetT->Fill(ztrackjTr);
09fce7b3 645
931f8b00 646 }//end loop on jet tracks
647
0b7f0a99 648 fhNtrkjzNok->Fill(zg1jtrk);
09fce7b3 649 Double_t diffpt=origPtJet-sumPtTracks;
650 Double_t diffpz=jet->Pz()-sumPzTracks;
651 if(zg1jtrk>0){
0b7f0a99 652 fhDiffPtTrPtJzNok->Fill(diffpt);
653 fhDiffPzTrPtJzNok->Fill(diffpz);
09fce7b3 654
655 }else{
0b7f0a99 656 fhDiffPtTrPtJzok->Fill(diffpt);
657 fhDiffPzTrPtJzok->Fill(diffpz);
09fce7b3 658 }
659
d9126a4c 660 if(candidates==0){
5a51a395 661
0b7f0a99 662 fhPtJetPerEvNoD->Fill(fPtJet);
d9126a4c 663 }
0455d618 664 if(!fJetOnlyMode) {
b2705b43 665 //Printf("N candidates %d ", candidates);
666 for(Int_t ic = 0; ic < candidates; ic++) {
0b7f0a99 667 fhstat->Fill(7);
b2705b43 668 // D* candidates
669 AliVParticle* charm=0x0;
670 charm=(AliVParticle*)fCandidateArray->At(ic);
671 if(!charm) continue;
672 AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
673 fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
674 if (fIsDInJet) FlagFlavour(jet);
0b7f0a99 675 if (jet->TestFlavourTag(AliEmcalJet::kDStar) || jet->TestFlavourTag(AliEmcalJet::kD0)) fhstat->Fill(4);
b2705b43 676
677 //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.
678
679 Double_t pjet[3];
680 jet->PxPyPz(pjet);
d8324fea 681 //It corrects the jet momentum if it was asked for jet background subtraction
60e4f65e 682 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
683 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
684 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 685
686 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 687 RecalculateMomentum(pjet,fPmissing);
d8324fea 688 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
b2705b43 689
09fce7b3 690
691 //debugging histograms
d8324fea 692 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
0b7f0a99 693 for(Int_t k=0;k<3;k++) fhMissPi[k]->Fill(fPmissing[k]);
b2705b43 694
0b7f0a99 695 fhmissingp->Fill(pmissing);
b2705b43 696 Double_t ptdiff=fPtJet-origPtJet;
0b7f0a99 697 fhDeltaPtJet->Fill(ptdiff);
698 fhRelDeltaPtJet->Fill(ptdiff/origPtJet);
b2705b43 699
700 FillHistogramsRecoJetCorr(charm, jet, aodEvent);
701
702 }//end loop on candidates
931f8b00 703
b2705b43 704 //retrieve side band background candidates for Dstar
705 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
931f8b00 706
b2705b43 707 for(Int_t ib=0;ib<nsbcand;ib++){
708 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
709 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
710 if(!sbcand) continue;
ce11164d 711
b2705b43 712 fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
713 Double_t pjet[3];
714 jet->PxPyPz(pjet);
d8324fea 715 //It corrects the jet momentum if it was asked for jet background subtraction
60e4f65e 716 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
717 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
718 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 719
720 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 721 RecalculateMomentum(pjet,fPmissing);
d8324fea 722 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
b2705b43 723
724 SideBandBackground(sbcand,jet);
725
726 }
727 if(fUseMCInfo){
da94b30f 728
b2705b43 729 AliAODRecoDecayHF* charmbg = 0x0;
a1a0a01b 730 charmbg=(AliAODRecoDecayHF*)fSideBandArray->At(ib);
b2705b43 731 if(!charmbg) continue;
0b7f0a99 732 fhstat->Fill(8);
b2705b43 733 fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
da94b30f 734 if (fIsDInJet) {
735 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 736 fhstat->Fill(9);
da94b30f 737 }
b2705b43 738 Double_t pjet[3];
739 jet->PxPyPz(pjet);
d8324fea 740 //It corrects the jet momentum if it was asked for jet background subtraction
a1a0a01b 741 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
742 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
743 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
d8324fea 744
745 //It corrects the jet momentum due to D daughters out of the jet
b2705b43 746 RecalculateMomentum(pjet,fPmissing);
d8324fea 747 fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
b2705b43 748
a1a0a01b 749 MCBackground(charmbg,jet);
b2705b43 750 }
931f8b00 751 }
752 }
753 } // end of jet loop
754
0b7f0a99 755 fhNJetPerEv->Fill(cntjet);
756 if(candidates==0) fhNJetPerEvNoD->Fill(cntjet);
76bf81f2 757 PostData(1,fOutput);
b5d0ba13 758 return kTRUE;
931f8b00 759}
760
761//_______________________________________________________________________________
762
763void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
764{
765 // The Terminate() function is the last function to be called during
766 // a query. It always runs on the client, it can be used to present
767 // the results graphically or save the results to file.
768
769 Info("Terminate"," terminate");
770 AliAnalysisTaskSE::Terminate();
771
76bf81f2 772 fOutput = dynamic_cast<TList*> (GetOutputData(1));
773 if (!fOutput) {
774 printf("ERROR: fOutput not available\n");
931f8b00 775 return;
776 }
777}
778
779//_______________________________________________________________________________
780
781void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
782 Float_t mass=0;
783 Int_t abspdg=TMath::Abs(pdg);
784
785 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
786 // compute the Delta mass for the D*
787 if(fCandidateType==kDstartoKpipi){
788 Float_t mass1=0;
789 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
790 mass = mass-mass1;
791 }
792
793 fMinMass = mass-range;
794 fMaxMass = mass+range;
795
796 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
797 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
798}
799
800//_______________________________________________________________________________
801
802void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
803 if(uplimit>lowlimit)
804 {
805 fMinMass = lowlimit;
806 fMaxMass = uplimit;
807 }
808 else{
809 printf("Error! Lower limit larger than upper limit!\n");
810 fMinMass = uplimit - uplimit*0.2;
811 fMaxMass = uplimit;
812 }
813}
814
815//_______________________________________________________________________________
816
817Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
818 if(nptbins>30) {
819 AliInfo("Maximum number of bins allowed is 30!");
820 return kFALSE;
821 }
822 if(!width) return kFALSE;
823 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
824 return kTRUE;
825}
826
827//_______________________________________________________________________________
828
09fce7b3 829Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
931f8b00 830 if(!part || !jet){
831 printf("Particle or jet do not exist!\n");
832 return -99;
833 }
834 Double_t p[3],pj[3];
835 Bool_t okpp=part->PxPyPz(p);
836 Bool_t okpj=jet->PxPyPz(pj);
60e4f65e 837
838 //Background Subtraction
d8324fea 839 //It corrects the each component of the jet momentum for Z calculation
60e4f65e 840 AliRhoParameter *rho = 0;
841 Double_t rhoval = 0;
842 TString sname("Rho");
843 if (!sname.IsNull()) {
844 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
845 if(rho){
846 rhoval = rho->GetVal();
847 pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
848 pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
849 pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
850 }
851 }
852
853
854
931f8b00 855 if(!okpp || !okpj){
856 printf("Problems getting momenta\n");
857 return -999;
858 }
b2705b43 859
860 RecalculateMomentum(pj, fPmissing);
09fce7b3 861 if(transverse) return ZT(p,pj);
862 else
863 return Z(p,pj);
864}
865
866//_______________________________________________________________________________
867Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
868
b2705b43 869 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
870 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
931f8b00 871 return z;
872}
873
09fce7b3 874
875//_______________________________________________________________________________
876Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
877
878 Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
879 Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
880 return z;
881}
882
931f8b00 883//_______________________________________________________________________________
884
b2705b43 885void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
886
887 pj[0]+=pmissing[0];
888 pj[1]+=pmissing[1];
889 pj[2]+=pmissing[2];
890
891}
892
893//_______________________________________________________________________________
894
931f8b00 895Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
896
897 // Statistics
da94b30f 898 Int_t nbins=8;
899 if(fUseMCInfo) nbins+=2;
900
0b7f0a99 901 fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
902 fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
903 fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
904 fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
905 fhstat->GetXaxis()->SetBinLabel(4,"N jets");
906 fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
907 fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
908 fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
909 fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
da94b30f 910 if(fUseMCInfo) {
0b7f0a99 911 fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
912 fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
913 fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
914 fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
da94b30f 915 }
0b7f0a99 916 fhstat->SetNdivisions(1);
917 fOutput->Add(fhstat);
931f8b00 918
e278882f 919 const Int_t nbinsmass=300;
931f8b00 920 const Int_t nbinsptjet=500;
921 const Int_t nbinsptD=100;
922 const Int_t nbinsz=100;
923 const Int_t nbinsphi=200;
1353fa73 924 const Int_t nbinseta=100;
ad6abcae 925
926 //binning for THnSparse
927 const Int_t nbinsSpsmass=50;
928 const Int_t nbinsSpsptjet=100;
929 const Int_t nbinsSpsptD=50;
930 const Int_t nbinsSpsz=100;
931 const Int_t nbinsSpsphi=100;
932 const Int_t nbinsSpseta=60;
933 const Int_t nbinsSpsContrib=100;
934 const Int_t nbinsSpsA=100;
935
931f8b00 936 const Float_t ptjetlims[2]={0.,200.};
937 const Float_t ptDlims[2]={0.,50.};
938 const Float_t zlims[2]={0.,1.2};
939 const Float_t philims[2]={0.,6.3};
1353fa73 940 const Float_t etalims[2]={-1.5,1.5};
0455d618 941 const Int_t nContriblims[2]={0,100};
942 const Float_t arealims[2]={0.,2};
1353fa73 943
931f8b00 944 // jet related fistograms
945
0b7f0a99 946 fhEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
947 fhEjetTrks->Sumw2();
948 fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
949 fhPhiJetTrks->Sumw2();
950 fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
951 fhEtaJetTrks->Sumw2();
952 fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
953 fhPtJetTrks->Sumw2();
954
955 fhEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
956 fhEjet->Sumw2();
957 fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
958 fhPhiJet->Sumw2();
959 fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
960 fhEtaJet->Sumw2();
961 fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
962 fhPtJet->Sumw2();
963
964 fhdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
965 fhdeltaRJetTracks->Sumw2();
966
967 fhNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
968 fhNtrArr->Sumw2();
969
970 fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
971 fhNJetPerEv->Sumw2();
972
e278882f 973 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
974 //NB at the moment fillet with jet track imp par!!!
975 fhVtxX = new TH1F("hVtxX", "Vertex X;vtx_x",100,-0.5,0.5);
976 fhVtxY = new TH1F("hVtxY", "Vertex Y;vtx_y",100,-0.5,0.5);
977 fhVtxZ = new TH1F("hVtxZ", "Vertex Z;vtx_z",100,-10,10);
978 /*
979 if(fUseMCInfo){
980 //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
981 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
982 fOutput->Add(fhImpParB);
983
984 }
985 */
0b7f0a99 986
987 fOutput->Add(fhEjetTrks);
988 fOutput->Add(fhPhiJetTrks);
989 fOutput->Add(fhEtaJetTrks);
990 fOutput->Add(fhPtJetTrks);
991 fOutput->Add(fhEjet);
992 fOutput->Add(fhPhiJet);
993 fOutput->Add(fhEtaJet);
994 fOutput->Add(fhPtJet);
995 fOutput->Add(fhdeltaRJetTracks);
996 fOutput->Add(fhNtrArr);
997 fOutput->Add(fhNJetPerEv);
e278882f 998 fOutput->Add(fhImpPar);
999 fOutput->Add(fhVtxX);
1000 fOutput->Add(fhVtxY);
1001 fOutput->Add(fhVtxZ);
ad6abcae 1002 if(fJetOnlyMode && fSwitchOnSparses){
0455d618 1003 //thnsparse for jets
1004 const Int_t nAxis=6;
ad6abcae 1005 const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
2942f542 1006 const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
0455d618 1007 const Double_t
2942f542 1008 maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
0b7f0a99 1009 fhsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
1010 fhsJet->Sumw2();
931f8b00 1011
0b7f0a99 1012 fOutput->Add(fhsJet);
0455d618 1013
931f8b00 1014 }
931f8b00 1015
0455d618 1016 if(!fJetOnlyMode){
09fce7b3 1017
1018 //debugging histograms
0b7f0a99 1019 fhControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
1020 fhControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
1021 fhControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
1022 fhControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
1023 fhControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
1024 fhControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
1025 fhControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
1026 fhControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
1027 fhControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
09fce7b3 1028
0b7f0a99 1029 fhControlDInJ->SetNdivisions(1);
1030 fhControlDInJ->GetXaxis()->SetLabelSize(0.05);
1031 fOutput->Add(fhControlDInJ);
09fce7b3 1032
0b7f0a99 1033 fhmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
1034 fOutput->Add(fhmissingp);
09fce7b3 1035
0b7f0a99 1036 fhMissPi=new TH1F*[3];
09fce7b3 1037 for(Int_t k=0;k<3;k++) {
0b7f0a99 1038 fhMissPi[k]=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
1039 fOutput->Add(fhMissPi[k]);
09fce7b3 1040 }
0b7f0a99 1041 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 1042
0b7f0a99 1043 fOutput->Add(fhDeltaPtJet);
1044 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.);
1045 fOutput->Add(fhRelDeltaPtJet);
09fce7b3 1046
0b7f0a99 1047 fhzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
1048 fOutput->Add(fhzDinjet);
09fce7b3 1049 //frag func of tracks in the jet
0b7f0a99 1050 fhztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
1051 fOutput->Add(fhztracksinjet);
09fce7b3 1052
1053 //check jet and tracks
0b7f0a99 1054 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);
1055 fOutput->Add(fhDiffPtTrPtJzok);
1056 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);
1057 fOutput->Add(fhDiffPtTrPtJzNok);
1058 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);
1059 fOutput->Add(fhDiffPzTrPtJzok);
1060 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);
1061 fOutput->Add(fhDiffPzTrPtJzNok);
1062 fhNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
1063 fOutput->Add(fhNtrkjzNok);
09fce7b3 1064
1065 //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
0b7f0a99 1066 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]);
1067 fOutput->Add(fhzDT);
1068 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]);
1069 fOutput->Add(fhztracksinjetT);
09fce7b3 1070
0b7f0a99 1071 fhIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
1072 fOutput->Add(fhIDddaugh);
1073 fhIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
1074 fOutput->Add(fhIDddaughOut);
1075 fhIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
1076 fOutput->Add(fhIDjetTracks);
09fce7b3 1077
0b7f0a99 1078 fhDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
1079 fOutput->Add(fhDRdaughOut);
09fce7b3 1080
1081
0455d618 1082 if(fCandidateType==kDstartoKpipi)
1083 {
da94b30f 1084 if(fSwitchOnSB){
0b7f0a99 1085 fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1086 fhDiffSideBand->SetStats(kTRUE);
1087 fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
1088 fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1089 fhDiffSideBand->Sumw2();
1090 fOutput->Add(fhDiffSideBand);
da94b30f 1091 }
0455d618 1092
0b7f0a99 1093 fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
1094 fhPtPion->SetStats(kTRUE);
1095 fhPtPion->GetXaxis()->SetTitle("GeV/c");
1096 fhPtPion->GetYaxis()->SetTitle("Entries");
1097 fhPtPion->Sumw2();
1098 fOutput->Add(fhPtPion);
0455d618 1099
1100 }
1101 // D related histograms
0b7f0a99 1102 fhNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
1103 fhNDPerEvNoJet->Sumw2();
1104 fOutput->Add(fhNDPerEvNoJet);
0455d618 1105
0b7f0a99 1106 fhptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
1107 fhptDPerEvNoJet->Sumw2();
1108 fOutput->Add(fhptDPerEvNoJet);
0455d618 1109
ad6abcae 1110 if(fSwitchOnSparses){
1111 const Int_t nAxisD=4;
1112 const Int_t nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
1113 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
1114 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
0b7f0a99 1115 fhsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
1116 fhsDstandalone->Sumw2();
ad6abcae 1117
0b7f0a99 1118 fOutput->Add(fhsDstandalone);
ad6abcae 1119 }
0455d618 1120
ad6abcae 1121 /*
0455d618 1122 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]);
1123 hPtJetWithD->Sumw2();
1124 //for the MC this histogram is filled with the real background
1125 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]);
1126 hPtJetWithDsb->Sumw2();
1127
ad6abcae 1128 fOutput->Add(hPtJetWithD);
1129 fOutput->Add(hPtJetWithDsb);
1130
1131 */
0b7f0a99 1132 fhNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
1133 fhNJetPerEvNoD->Sumw2();
0455d618 1134
0b7f0a99 1135 fhPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1136 fhPtJetPerEvNoD->Sumw2();
0455d618 1137
0b7f0a99 1138 fOutput->Add(fhNJetPerEvNoD);
1139 fOutput->Add(fhPtJetPerEvNoD);
0455d618 1140
0b7f0a99 1141 fhDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1142 fhDeltaRD->Sumw2();
1143 fOutput->Add(fhDeltaRD);
0455d618 1144
0b7f0a99 1145 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]);
1146 fhDeltaRptDptj->Sumw2();
1147 fOutput->Add(fhDeltaRptDptj);
5a51a395 1148
1149 if(fUseMCInfo){
0b7f0a99 1150 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]);
1151 fhDeltaRptDptjB->Sumw2();
1152 fOutput->Add(fhDeltaRptDptjB);
5a51a395 1153 }
1154
0455d618 1155 //background (side bands for the Dstar and like sign for D0)
cc5dbb3b 1156 AliJetContainer *jetCont=GetJetContainer(0);
1157 if(!jetCont){
1158 Printf("Container 0 not found, try with name %s", fJetArrName.Data());
1159 jetCont=GetJetContainer(fJetArrName);
1160 if(!jetCont) Printf("NOT FOUND AGAIN");
1161 return kFALSE;
1162 }
1163 Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
1164 //fJetRadius=jetCont->GetJetRadius();
0b7f0a99 1165 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]);
1166 fhInvMassptD->SetStats(kTRUE);
1167 fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1168 fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1169 fhInvMassptD->Sumw2();
0455d618 1170
0b7f0a99 1171 fOutput->Add(fhInvMassptD);
0455d618 1172
1173 if(fUseMCInfo){
0b7f0a99 1174 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]);
1175 fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1176 fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1177 fhInvMassptDbg->Sumw2();
1178 fOutput->Add(fhInvMassptDbg);
0455d618 1179
1180 }
0455d618 1181
ad6abcae 1182 if(fSwitchOnSparses){
1183 Double_t pi=TMath::Pi(), philims2[2];
1184 philims2[0]=-pi*1./2.;
1185 philims2[1]=pi*3./2.;
0b7f0a99 1186 fhsDphiz=0x0; //definition below according to the switches
ad6abcae 1187
1188 if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1189 AliInfo("Creating a 9 axes container: This might cause large memory usage");
1190 const Int_t nAxis=9;
1191 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1192 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1193 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1194 fNAxesBigSparse=nAxis;
1195
0b7f0a99 1196 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 1197 }
1198
1199 if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1200 fSwitchOnPhiAxis=0;
1201 fSwitchOnOutOfConeAxis=0;
1202 fSwitchOnSB=0;
a1a0a01b 1203 if(fUseMCInfo){
1204 AliInfo("Creating a 7 axes container (MB background candidates)");
da94b30f 1205 const Int_t nAxis=7;
a1a0a01b 1206 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1207 const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1208 const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1209 fNAxesBigSparse=nAxis;
0b7f0a99 1210 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 1211
1212 } else{
1213 AliInfo("Creating a 6 axes container");
1214 const Int_t nAxis=6;
1215 const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1216 const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1217 const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1218 fNAxesBigSparse=nAxis;
1219
0b7f0a99 1220 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 1221 }
ad6abcae 1222 }
0b7f0a99 1223 if(!fhsDphiz) AliFatal("No THnSparse created");
1224 fhsDphiz->Sumw2();
ad6abcae 1225
0b7f0a99 1226 fOutput->Add(fhsDphiz);
ad6abcae 1227 }
0455d618 1228 }
76bf81f2 1229 PostData(1, fOutput);
931f8b00 1230
1231 return kTRUE;
1232}
1233
1234//_______________________________________________________________________________
1235
8577f7bc 1236void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet, AliAODEvent* aodEvent){
931f8b00 1237
1238 Double_t ptD=candidate->Pt();
d8324fea 1239 Double_t deltaR=DeltaR(jet,candidate);
b2705b43 1240 Double_t phiD=candidate->Phi();
1241 Double_t deltaphi = jet->Phi()-phiD;
1242 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1243 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
931f8b00 1244 Double_t z=Z(candidate,jet);
09fce7b3 1245 /*
1246 if(z>1) {
1247 ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1248 Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1249
1250 for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1251
1252 ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1253 Double_t ptdiff=fPtJet-jet->Pt();
1254 ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1255 ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1256
1257
1258 }
1259 */
0b7f0a99 1260 if(fIsDInJet)fhzDT->Fill(Z(candidate,jet,kTRUE));
1261
1262 fhDeltaRD->Fill(deltaR);
1263 fhDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
5a51a395 1264
ce11164d 1265 Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1266 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
931f8b00 1267 if(fUseReco){
1268 if(fCandidateType==kD0toKpi) {
1269 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
8577f7bc 1270
ad6abcae 1271 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
931f8b00 1272
1273 }
1274
1275 if(fCandidateType==kDstartoKpipi) {
1276 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
ad6abcae 1277 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
931f8b00 1278
1279 }
1280 } else{ //generated level
1281 //AliInfo("Non reco");
ad6abcae 1282 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
931f8b00 1283 }
1284}
1285
1286//_______________________________________________________________________________
1287
ad6abcae 1288void 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 1289
931f8b00 1290
1291 Double_t masses[2]={0.,0.};
1292 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1293 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1294
1295 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1296 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1297
ad6abcae 1298 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
ad6abcae 1299 Double_t *point=0x0;
1300 if(fNAxesBigSparse==9){
1301 point=new Double_t[9];
1302 point[0]=z;
1303 point[1]=dPhi;
1304 point[2]=ptj;
1305 point[3]=ptD;
1306 point[4]=masses[0];
1307 point[5]=0;
1308 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1309 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1310 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1311 }
1312 if(fNAxesBigSparse==6){
1313 point=new Double_t[6];
1314 point[0]=z;
1315 point[1]=ptj;
1316 point[2]=ptD;
1317 point[3]=masses[0];
1318 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1319 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
a1a0a01b 1320}
1321 if(fNAxesBigSparse==7){
1322 point=new Double_t[7];
1323 point[0]=z;
1324 point[1]=ptj;
1325 point[2]=ptD;
1326 point[3]=masses[0];
1327 point[4]=0;
1328 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1329 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1330
ad6abcae 1331 }
1332
1333
e278882f 1334 //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
931f8b00 1335 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1336 if(isselected==1 || isselected==3) {
1337
ad6abcae 1338 //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
931f8b00 1339
b2705b43 1340 FillMassHistograms(masses[0], ptD);
0b7f0a99 1341 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1342 }
1343 if(isselected>=2) {
ad6abcae 1344 //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
931f8b00 1345
b2705b43 1346 FillMassHistograms(masses[1], ptD);
ad6abcae 1347 if(fNAxesBigSparse==9) point[4]=masses[1];
a1a0a01b 1348 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
0b7f0a99 1349 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1350 }
1351
1352}
1353
1354//_______________________________________________________________________________
1355
ad6abcae 1356void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
931f8b00 1357 //dPhi and z not used at the moment,but will be (re)added
1358
1359 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1360 Double_t deltamass= dstar->DeltaInvMass();
1361 //Double_t massD0= dstar->InvMassD0();
1362
1363
0b7f0a99 1364 fhPtPion->Fill(softpi->Pt());
931f8b00 1365
ad6abcae 1366 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
ed5a88ff 1367 Int_t isSB=0;//IsDzeroSideBand(dstar);
931f8b00 1368
ad6abcae 1369 //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1370 Double_t *point=0x0;
1371 if(fNAxesBigSparse==9){
1372 point=new Double_t[9];
1373 point[0]=z;
1374 point[1]=dPhi;
1375 point[2]=ptj;
1376 point[3]=ptD;
1377 point[4]=deltamass;
1378 point[5]=static_cast<Double_t>(isSB);
1379 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1380 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1381 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1382 }
1383 if(fNAxesBigSparse==6){
1384 point=new Double_t[6];
1385 point[0]=z;
1386 point[1]=ptj;
1387 point[2]=ptD;
1388 point[3]=deltamass;
1389 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1390 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1391 }
a1a0a01b 1392 if(fNAxesBigSparse==7){
1393 point=new Double_t[7];
1394 point[0]=z;
1395 point[1]=ptj;
1396 point[2]=ptD;
1397 point[3]=deltamass;
1398 point[4]=0;
1399 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1400 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1401 }
931f8b00 1402
ad6abcae 1403 //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
931f8b00 1404
b2705b43 1405 FillMassHistograms(deltamass, ptD);
0b7f0a99 1406 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1407
1408}
1409
1410//_______________________________________________________________________________
1411
ad6abcae 1412void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
931f8b00 1413
1414 Double_t pdgmass=0;
a1a0a01b 1415 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1416 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
ad6abcae 1417 //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
0b7f0a99 1418
ad6abcae 1419 //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
5a51a395 1420
ad6abcae 1421 Double_t *point=0x0;
1422 if(fNAxesBigSparse==9){
1423 point=new Double_t[9];
1424 point[0]=z;
1425 point[1]=dPhi;
1426 point[2]=ptjet;
1427 point[3]=ptD;
1428 point[4]=pdgmass;
1429 point[5]=0;
1430 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1431 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1432 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1433 }
1434 if(fNAxesBigSparse==6){
1435 point=new Double_t[6];
1436 point[0]=z;
1437 point[1]=ptjet;
1438 point[2]=ptD;
1439 point[3]=pdgmass;
1440 point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1441 point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1442 }
a1a0a01b 1443 if(fNAxesBigSparse==7){
da94b30f 1444 point=new Double_t[7];
a1a0a01b 1445 point[0]=z;
1446 point[1]=ptjet;
1447 point[2]=ptD;
1448 point[3]=pdgmass;
1449 point[4]=1;
1450 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1451 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1452 }
931f8b00 1453
a1a0a01b 1454
1455
ad6abcae 1456 if(fNAxesBigSparse==9) point[4]=pdgmass;
a1a0a01b 1457 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
0b7f0a99 1458 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
ad6abcae 1459 //if(fIsDInJet) {
1460 // hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1461 //}
931f8b00 1462
1463}
1464
1465//_______________________________________________________________________________
1466
b2705b43 1467void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
931f8b00 1468
b2705b43 1469 if(fIsDInJet) {
0b7f0a99 1470 fhInvMassptD->Fill(mass,ptD);
931f8b00 1471 }
1472}
1473
1353fa73 1474//________________________________________________________________________________
1475
b2705b43 1476void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1477
931f8b00 1478 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1479 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
b2705b43 1480 if (fIsDInJet) jet->AddFlavourTag(tag);
931f8b00 1481
1482 return;
1483
1484}
1485
1486//_______________________________________________________________________________
1487
1488void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1489
1490 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1491 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1492 // has a width of ~5 sigmas. Two band needed for opening angle considerations
0b7f0a99 1493
ad6abcae 1494 //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
931f8b00 1495
ce11164d 1496 Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1497 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1498
931f8b00 1499 Double_t deltaM=candDstar->DeltaInvMass();
1500 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
ed5a88ff 1501 Double_t z=Z(candDstar,jet);
931f8b00 1502 Double_t ptD=candDstar->Pt();
b2705b43 1503
931f8b00 1504 Double_t dPhi=jet->Phi()-candDstar->Phi();
b2705b43 1505
1506 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1507 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
931f8b00 1508
ed5a88ff 1509 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
ad6abcae 1510 //if no SB analysis we should not enter here, so no need of checking the number of axes
1511 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 1512 //fill the background histogram with the side bands or when looking at MC with the real background
1513 if(isSideBand==1){
0b7f0a99 1514 fhDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
931f8b00 1515 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
0b7f0a99 1516 if(fSwitchOnSparses) fhsDphiz->Fill(point,1.);
ad6abcae 1517 //if(fIsDInJet){ // evaluate in the near side
1518 // hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1519 //}
931f8b00 1520 }
1521}
1522
1523//_______________________________________________________________________________
1524
a1a0a01b 1525void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
931f8b00 1526
1527 //need mass, deltaR, pt jet, ptD
1528 //two cases: D0 and Dstar
1529
1530 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1531 if(!isselected) return;
1532
1533 Double_t ptD=candbg->Pt();
d8324fea 1534 Double_t deltaR=DeltaR(jet,candbg);
a1a0a01b 1535 Double_t phiD=candbg->Phi();
1536 Double_t deltaphi = jet->Phi()-phiD;
1537 if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1538 if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1539 Double_t z=Z(candbg,jet);
1540
0b7f0a99 1541 if(fIsDInJet) fhzDT->Fill(Z(candbg,jet,kTRUE));
da94b30f 1542
5a51a395 1543
0b7f0a99 1544
1545
1546 fhDeltaRD->Fill(deltaR);
1547 fhDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
da94b30f 1548
a1a0a01b 1549 Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1550 Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1551
ad6abcae 1552 //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
ce11164d 1553
a1a0a01b 1554 Double_t *point=0x0;
1555 if(fNAxesBigSparse==9){
1556 point=new Double_t[9];
1557 point[0]=z;
1558 point[1]=deltaphi;
1559 point[2]=fPtJet;
1560 point[3]=ptD;
1561 point[4]=-999; //set below
1562 point[5]=1;
1563 point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1564 point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1565 point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1566 }
1567
1568 if(fNAxesBigSparse==7){
1569 point=new Double_t[7];
1570 point[0]=z;
1571 point[1]=fPtJet;
1572 point[2]=ptD;
1573 point[3]=-999; //set below
1574 point[4]=1;
1575 point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1576 point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1577 }
ce11164d 1578
931f8b00 1579 if(fCandidateType==kDstartoKpipi){
1580 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1581 Double_t deltaM=dstarbg->DeltaInvMass();
0b7f0a99 1582 fhInvMassptDbg->Fill(deltaM,ptD);
ad6abcae 1583 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
a1a0a01b 1584 if(fNAxesBigSparse==9) point[4]=deltaM;
1585 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
0b7f0a99 1586 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
931f8b00 1587 }
1588
1589 if(fCandidateType==kD0toKpi){
1590 Double_t masses[2]={0.,0.};
1591 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1592 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1593
1594 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1595 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1596
1597
1598 if(isselected==1 || isselected==3) {
ad6abcae 1599 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
0b7f0a99 1600 fhInvMassptDbg->Fill(masses[0],ptD);
a1a0a01b 1601 if(fNAxesBigSparse==9) point[4]=masses[0];
1602 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
0b7f0a99 1603 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
a1a0a01b 1604 }
931f8b00 1605 if(isselected>=2) {
ad6abcae 1606 //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
0b7f0a99 1607 fhInvMassptDbg->Fill(masses[1],ptD);
a1a0a01b 1608 if(fNAxesBigSparse==9) point[4]=masses[1];
1609 if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
0b7f0a99 1610 if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
a1a0a01b 1611
931f8b00 1612 }
1613
1614
1615 }
1616}
1617
1618//_______________________________________________________________________________
1619
d8324fea 1620Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
931f8b00 1621 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
d8324fea 1622 //It recalculates the eta-phi values if it was asked for background subtraction of the jets
931f8b00 1623 if(!p1 || !p2) return -1;
d8324fea 1624
1625 Double_t phi1=p1->Phi(),eta1=p1->Eta();
1626
1627 //It subtracts the backgroud of jets if it was asked for it.
1628 TString sname("Rho");
1629 if (!sname.IsNull()) {
1630 AliRhoParameter *rho = 0;
1631 rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
1632 if(rho){
1633 Double_t pj[3];
1634 Bool_t okpj=p1->PxPyPz(pj);
1635 if(!okpj){
1636 printf("Problems getting momenta\n");
1637 return -999;
1638 }
1639 Double_t rhoval = rho->GetVal();
1640 pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1641 pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1642 pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1643 //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1644 //[0,pi] if py > 0 and
1645 //[pi,2pi] if py < 0
1646 phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1647 if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1648 eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1649 }
1650 }
1651
931f8b00 1652 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1653
1654 Double_t dPhi=phi1-phi2;
1655 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1656 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1657
1658 Double_t dEta=eta1-eta2;
1659 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1660 return deltaR;
1661
1662}
1663
1664//_______________________________________________________________________________
1665
1666Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1667
1668 Double_t ptD=candDstar->Pt();
1669 Int_t bin = fCuts->PtBin(ptD);
1670 if (bin < 0){
1671 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1672 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1673 return -1; // out of bounds
1674 }
1675
1676 Double_t invM=candDstar->InvMassD0();
1677 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1678
1679 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1680 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1681
1682 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1683 else return 0;
1684
1685}
b2705b43 1686
1687//_______________________________________________________________________________
1688
1689Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1690 //returns true/false and the array of particles not found among jet constituents
1691
b2705b43 1692
1693 Int_t daughtersID[3];
1694 Int_t ndaugh=0;
1695 Int_t dmatchedID[3]={0,0,0};
1696 Int_t countmatches=0;
1697 if (fCandidateType==kDstartoKpipi){
1698 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1699 AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1700 daughtersID[0]=dzero->GetProngID(0);
1701 daughtersID[1]=dzero->GetProngID(1);
1702 daughtersID[2]=dstar->GetBachelor()->GetID();
1703 ndaugh=3;
1704
1705 dtrks=new AliAODTrack*[3];
1706 dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1707 dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1708 dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1709
1710 //check
1711 if(fillH){
0b7f0a99 1712 if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID()) fhControlDInJ->Fill(6);
b2705b43 1713
0b7f0a99 1714 fhIDddaugh->Fill(daughtersID[0]);
1715 fhIDddaugh->Fill(daughtersID[1]);
1716 fhIDddaugh->Fill(daughtersID[2]);
b2705b43 1717
1718 }
1719 //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1720 }
1721
1722 if (fCandidateType==kD0toKpi){
1723 daughtersID[0]=charm->GetProngID(0);
1724 daughtersID[1]=charm->GetProngID(1);
1725 ndaugh=2;
1726 if(fillH){
0b7f0a99 1727 fhIDddaugh->Fill(daughtersID[0]);
1728 fhIDddaugh->Fill(daughtersID[1]);
b2705b43 1729 }
1730 dtrks=new AliAODTrack*[2];
1731 dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1732 dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1733
1734 }
1735
1736 const Int_t cndaugh=ndaugh;
1737 daughOutOfJetID=new Int_t[cndaugh];
1738
1739 Int_t nchtrk=thejet->GetNumberOfTracks();
1740 for(Int_t itk=0;itk<nchtrk;itk++){
1741 AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
e278882f 1742 if(!tkjet) continue;
b2705b43 1743 Int_t idtkjet=tkjet->GetID();
0b7f0a99 1744 if(fillH) fhIDjetTracks->Fill(idtkjet);
b2705b43 1745 for(Int_t id=0;id<ndaugh;id++){
1746 if(idtkjet==daughtersID[id]) {
1747 dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1748 countmatches++;
1749
1750 }
1751
1752 if(countmatches==ndaugh) break;
1753 }
1754 }
1755 //reverse: include in the array the ID of daughters not matching
1756
1757 for(Int_t id=0;id<ndaugh;id++){
1758 if(dmatchedID[id]==0) {
1759 daughOutOfJetID[id]=daughtersID[id];
0b7f0a99 1760 if(fillH) fhIDddaughOut->Fill(daughOutOfJetID[id]);
b2705b43 1761 }
1762 else daughOutOfJetID[id]=0;
1763 }
1764 if(fillH){
0b7f0a99 1765 if((ndaugh-countmatches) == 1) fhControlDInJ->Fill(0);
1766 if((ndaugh-countmatches) == 2) fhControlDInJ->Fill(1);
1767 if((ndaugh-countmatches) == 3) fhControlDInJ->Fill(2);
b2705b43 1768 }
1769 if(ndaugh!=countmatches){
1770 return kFALSE;
1771 }
1772
1773 return kTRUE;
1774}
1775
1776//_______________________________________________________________________________
1777
1778Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1779
1780 //check the conditions type:
1781 //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1782 //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1783 //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
e712e4fb 1784 //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 1785
09fce7b3 1786 fPmissing[0]=0;
1787 fPmissing[1]=0;
1788 fPmissing[2]=0;
b2705b43 1789
1790 Bool_t testDeltaR=kFALSE;
1791 if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1792
1793 Int_t* daughOutOfJet;
1794 AliAODTrack** charmDaugh;
1795 Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
09fce7b3 1796 if(testDaugh && testDeltaR) {
1797 //Z of candidates with daughters in the jet
0b7f0a99 1798 fhzDinjet->Fill(Z(charm,thejet));
09fce7b3 1799
1800 }
e712e4fb 1801 if(!testDaugh && testDeltaR && fTypeDInJet>=2){
b2705b43 1802
1803 Int_t ndaugh=3;
1804 if(fCandidateType==kD0toKpi) ndaugh=2;
1805 Int_t nOut=ndaugh;
1806
b2705b43 1807 for(Int_t id=0;id<ndaugh;id++){
1808 if(daughOutOfJet[id]!=0){ //non-matched daughter
1809 //get track and its momentum
1810 nOut--;
1811 if(fillH) {
0b7f0a99 1812 fhControlDInJ->Fill(3);
1813 if(id<2) fhControlDInJ->Fill(4);
1814 if(id==2)fhControlDInJ->Fill(5);
1815 fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
b2705b43 1816 }
e712e4fb 1817 if(fTypeDInJet==2){
1818 fPmissing[0]+=charmDaugh[id]->Px();
1819 fPmissing[1]+=charmDaugh[id]->Py();
1820 fPmissing[2]+=charmDaugh[id]->Pz();
1821 }
1822 if(fTypeDInJet==3){
1823 Double_t ptdaug = charmDaugh[id]->Pt();
1824 Double_t ptjet = thejet->Pt();
1825 Double_t ptn = ptjet+ptdaug;
1826 Double_t phidaug = charmDaugh[id]->Phi();
1827 Double_t phijet = thejet->Phi();
1828 Double_t phin = (phijet/ptjet+phidaug/ptdaug)/(1./ptjet+ 1./ptdaug);
1829 Double_t etadaug = charmDaugh[id]->Eta();
1830 Double_t etajet = thejet->Eta();
1831 Double_t etan = (etajet/ptjet+etadaug/ptdaug)/(1./ptjet+ 1./ptdaug);
1832
1833 fPmissing[0]+= ptn*TMath::Cos(phin);
1834 fPmissing[1]+= ptn*TMath::Sin(phin);
1835 fPmissing[2]+= ptn*TMath::SinH(etan);
1836 }
b2705b43 1837 }
1838
1839 }
1840
1841 //now the D in within the jet
1842 testDaugh=kTRUE;
1843
1844 }
1845
1846 delete[] charmDaugh;
1847
1848 Bool_t result=0;
1849 switch(fTypeDInJet){
1850 case 0:
1851 result=testDeltaR;
09fce7b3 1852 break;
b2705b43 1853 case 1:
1854 result=testDeltaR && testDaugh;
09fce7b3 1855 break;
e712e4fb 1856 case 2: //this case defines fPmissing
1857 result=testDeltaR && testDaugh;
09fce7b3 1858 break;
1859 default:
1860 AliInfo("Selection type not specified, use 1");
1861 result=testDeltaR && testDaugh;
1862 break;
b2705b43 1863 }
1864 return result;
1865}
ce11164d 1866
1867Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1868 //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1869
1870 Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1871 Bool_t binEMCal=kTRUE;
1872 Double_t phi=vpart->Phi(), eta=vpart->Eta();
1873 if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1874 if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1875 return binEMCal;
1876
1877
1878}