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