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