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