]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourFilterAndJetCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix mem leak from TClonesArray
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskFlavourFilterAndJetCorrelations.cxx
CommitLineData
ae040693 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>
33
34#include "AliAnalysisTaskFlavourFilterAndJetCorrelations.h"
35#include "AliAODMCHeader.h"
36#include "AliAODHandler.h"
37#include "AliAnalysisManager.h"
38#include "AliLog.h"
39#include "AliEmcalJet.h"
40#include "AliJetContainer.h"
41#include "AliAODRecoDecay.h"
42#include "AliAODRecoCascadeHF.h"
43#include "AliAODRecoDecayHF2Prong.h"
44#include "AliESDtrack.h"
45#include "AliAODMCParticle.h"
46#include "AliPicoTrack.h"
47#include "AliRDHFCutsD0toKpi.h"
48#include "AliRDHFCutsDStartoKpipi.h"
49
50ClassImp(AliAnalysisTaskFlavourFilterAndJetCorrelations)
51
52
53//_______________________________________________________________________________
54
55AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations() :
56AliAnalysisTaskEmcalJet("",kFALSE),
57fUseMCInfo(kTRUE),
58fUseReco(kTRUE),
59fCandidateType(),
60fPDGmother(),
61fNProngs(),
62fPDGdaughters(),
63fBranchName(),
64fmyOutput(0),
963e733a 65fmyOutputF(0),
ae040693 66fCuts(0),
67fMinMass(),
68fMaxMass(),
69fJetArrName(0),
70fCandArrName(0),
71fLeadingJetOnly(kFALSE),
963e733a 72fJetRadius(0),
73fCandidateArray(0),
74fSideBandArray(0)
ae040693 75{
76 //
77 // Default ctor
78 //
79 //SetMakeGeneralHistograms(kTRUE);
80
81}
82
83//_______________________________________________________________________________
84
85AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
86AliAnalysisTaskEmcalJet(name,kFALSE),
87fUseMCInfo(kTRUE),
88fUseReco(kTRUE),
89fCandidateType(),
90fPDGmother(),
91fNProngs(),
92fPDGdaughters(),
93fBranchName(),
94fmyOutput(0),
95fmyOutputF(0),
96fCuts(0),
97fMinMass(),
98fMaxMass(),
99fJetArrName(0),
100fCandArrName(0),
101fLeadingJetOnly(kFALSE),
963e733a 102fJetRadius(0),
103fCandidateArray(0),
104fSideBandArray(0)
ae040693 105{
106 //
107 // Constructor. Initialization of Inputs and Outputs
108 //
109
110 Info("AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Constructor");
111 fCuts=cuts;
112 fCandidateType=candtype;
113 const Int_t nptbins=fCuts->GetNPtBins();
114 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};
115
116 switch(fCandidateType){
117 case 0:
118 fPDGmother=421;
119 fNProngs=2;
120 fPDGdaughters[0]=211;//pi
121 fPDGdaughters[1]=321;//K
122 fPDGdaughters[2]=0; //empty
123 fPDGdaughters[3]=0; //empty
124 fBranchName="D0toKpi";
125 fCandArrName="D0";
126 break;
127 case 1:
128 fPDGmother=413;
129 fNProngs=3;
130 fPDGdaughters[1]=211;//pi soft
131 fPDGdaughters[0]=421;//D0
132 fPDGdaughters[2]=211;//pi fromD0
133 fPDGdaughters[3]=321; // K from D0
134 fBranchName="Dstar";
135 fCandArrName="Dstar";
136
137 if(nptbins<=13){
138 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
139 } else {
140 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
141 }
142 break;
143 default:
144 printf("%d not accepted!!\n",fCandidateType);
145 break;
146 }
147
148 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
149 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
150
151 DefineOutput(1,TList::Class()); // histos
152 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
153 DefineOutput(3,TList::Class()); //histos filter
154
155}
156
157//_______________________________________________________________________________
158
159AliAnalysisTaskFlavourFilterAndJetCorrelations::~AliAnalysisTaskFlavourFilterAndJetCorrelations() {
160 //
161 // destructor
162 //
163
164 Info("~AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Destructor");
165
166 delete fmyOutput;
167 delete fmyOutputF;
168 delete fCuts;
169
170}
171
172//_______________________________________________________________________________
173
174void AliAnalysisTaskFlavourFilterAndJetCorrelations::Init(){
175 //
176 // Initialization
177 //
178
179 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
180 switch(fCandidateType){
181 case 0:
182 {
183 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
184 copyfCuts->SetName("AnalysisCutsDzero");
185 // Post the data
186 PostData(2,copyfCuts);
187 }
188 break;
189 case 1:
190 {
191 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
192 copyfCuts->SetName("AnalysisCutsDStar");
193 // Post the cuts
194 PostData(2,copyfCuts);
195 }
196 break;
197 default:
198 return;
199 }
200
201 return;
202}
203
204//_______________________________________________________________________________
205
206void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserCreateOutputObjects() {
963e733a 207
208 fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
209 fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen"));
210
211 if (fCandidateType==kDstartoKpipi){
212 fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!
213 fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen"));
214 }
215
ae040693 216 // output
217 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
218 fmyOutput = new TList();
219 fmyOutput->SetOwner();
220 fmyOutput->SetName("pippo");
221 // define histograms
222 DefineHistoForAnalysis();
223 PostData(1,fmyOutput);
224
225 fmyOutputF = new TList();
226 fmyOutputF->SetOwner();
227 fmyOutputF->SetName("pluto");
228 // define histograms
229 DefineHistoForFiltering();
230 PostData(3,fmyOutputF);
231
232 return;
233}
234
235//_______________________________________________________________________________
236
237void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserExec(Option_t *)
238{
239 // user exec
240 if (!fInitialized){
241 AliAnalysisTaskEmcalJet::ExecOnce();
242 }
243
244 // Load the event
245 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
246
247 TClonesArray *arrayDStartoD0pi=0;
248 TClonesArray *trackArr = 0;
249
250 if (!aodEvent && AODEvent() && IsStandardAOD()) {
251
252 // In case there is an AOD handler writing a standard AOD, use the AOD
253 // event in memory rather than the input (ESD) event.
254 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
255
256 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
257 // have to taken from the AOD event hold by the AliAODExtension
258 AliAODHandler* aodHandler = (AliAODHandler*)
259 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
260 if(aodHandler->GetExtensions()) {
261 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
262 AliAODEvent *aodFromExt = ext->GetAOD();
263 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
264 }
265 } else {
266 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
267 }
268
269 if (!arrayDStartoD0pi) {
270 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
271 // return;
272 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
273
274 TClonesArray* mcArray = 0x0;
275 if (fUseMCInfo) {
276 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
277 if (!mcArray) {
278 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
279 return;
280 }
281 }
282
283 //retrieve jets
284 trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
285 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
286 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
287 fJetRadius=GetJetContainer(0)->GetJetRadius();
288
289 if(!trackArr){
290 AliInfo(Form("Could not find the track array\n"));
291 return;
292 }
293
294
295 //FILTER D MESONS
963e733a 296
297 fCandidateArray->Clear();
298 fSideBandArray->Clear();
ae040693 299
300 //Histograms
301 TH1I* hstatF = (TH1I*)fmyOutputF->FindObject("hstatF");
302 TH1F* hnSBCandEv=(TH1F*)fmyOutputF->FindObject("hnSBCandEv");
303 TH1F* hnCandEv=(TH1F*)fmyOutputF->FindObject("hnCandEv");
304 TH2F* hInvMassptDF = (TH2F*)fmyOutputF->FindObject("hInvMassptDF");
305
306 TH1F* hPtPion=0x0;
307 if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fmyOutputF->FindObject("hPtPion");
308 hstatF->Fill(0);
309
310 // fix for temporary bug in ESDfilter
311 // the AODs with null vertex pointer didn't pass the PhysSel
312 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
313
314 //Event selection
315 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
316 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
317 if (!iseventselected) return;
318 hstatF->Fill(1);
319
320 const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
321 if(fUseReco) hstatF->Fill(2, nD);
322
323 //D* and D0 prongs needed to MatchToMC method
324 Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
325 Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
326
327 //D0 from D0 bar
328 Int_t pdgdaughtersD0[2] = { 211, 321 }; // pi,K
329 Int_t pdgdaughtersD0bar[2] = { 321, 211 }; // K,pi
330
331 Int_t iCand =0;
332 Int_t iSBCand=0;
333 Int_t isSelected = 0;
334 AliAODRecoDecayHF *charmCand = 0;
335 AliAODMCParticle *charmPart = 0;
336 Bool_t isMCBkg=kFALSE;
337
338 Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
339
340 Int_t mcLabel = -9999;
341 Int_t pdgMeson = 413;
342 if (fCandidateType==kD0toKpi) pdgMeson = 421;
343
344 for (Int_t icharm=0; icharm<nD; icharm++) { //loop over D candidates
345 charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
346 if (!charmCand) continue;
347
348
349 if (fUseMCInfo) { // Look in MC, try to simulate the z
350 if (fCandidateType==kDstartoKpipi) {
351 AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
352 mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
353 }
354
355 if (fCandidateType==kD0toKpi)
356 mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
357
358 if (mcLabel<=0) isMCBkg=kTRUE;
359 else hstatF->Fill(2);
360 if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
361 }
362
363 Double_t ptD = charmCand->Pt();
364
365 // region of interest + cuts
366 if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
367
368 if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
369 //select by track cuts the side band candidates (don't want mass cut)
370 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
371 if (!isSelected) continue;
372 //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
373 Int_t bin = fCuts->PtBin(ptD);
374 if(bin<0 || bin>=fCuts->GetNPtBins()) {
375 AliError(Form("Pt %.3f out of bounds",ptD));
376 continue;
377 }
378 AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
379 //if data and Dstar from D0 side band
380 if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/|| ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
381
382 new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);
383 iSBCand++;
384 }
385 }
386 //candidate selected by cuts and PID
387 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
388 if (!isSelected) continue;
389
390 //for data and MC signal fill fCandidateArray
391 if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
392 // for data or MC with the requirement fUseReco fill with candidates
393 if(fUseReco) {
394 new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
395 //Printf("Filling reco");
396 hstatF->Fill(3);
397 }
398 // for MC with requirement particle level fill with AliAODMCParticle
399 else if (fUseMCInfo) {
400 new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
401 //Printf("Filling gen");
402 hstatF->Fill(3);
403 }
404
405 iCand++;
406 }
407 //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
408 else if(fUseReco){
409 new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
410 iSBCand++;
411 }
412
413
414 Double_t masses[2];
415 if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
416
417 //softpion from D* decay
418 AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
419 AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();
420
421 // select D* in the D0 window.
422 // In the cut object window is loose to allow for side bands
423
424
425 // retrieve the corresponding pt bin for the candidate
426 // and set the expected D0 width (x3)
427 // static const Int_t n = fCuts->GetNPtBins();
428 Int_t bin = fCuts->PtBin(ptD);
429 if(bin<0 || bin>=fCuts->GetNPtBins()) {
430 AliError(Form("Pt %.3f out of bounds",ptD));
431 continue;
432 }
433
434 AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
435 //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
436 if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
437 masses[0] = temp->DeltaInvMass(); //D*
438 masses[1] = 0.; //dummy for D*
439
440 //D* delta mass
441 hInvMassptDF->Fill(masses[0], ptD); // 2 D slice for pt bins
442
443 // D* pt and soft pion pt for good candidates
444 Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
445 Double_t invmassDelta = temp->DeltaInvMass();
446 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
447 }
448 } //Dstar specific
449
450 if (fCandidateType==kD0toKpi) { //D0->Kpi
451
452 //needed quantities
453 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
454 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
455 hstatF->Fill(3);
456
457 // mass vs pt
458 if (isSelected==1 || isSelected==3) hInvMassptDF->Fill(masses[0],ptD);
459 if (isSelected>=2) hInvMassptDF->Fill(masses[1],ptD);
460 } //D0 specific
461
462 charmCand = 0;
463 isMCBkg=kFALSE;
464 } // end of D cand loop
465
466 hnCandEv->Fill(fCandidateArray->GetEntriesFast());
467 if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
468 Int_t nsbcand=fSideBandArray->GetEntriesFast();
469 hstatF->Fill(4,nsbcand);
470 hnSBCandEv->Fill(nsbcand);
471 }
472
473 //CORRELATION WITH JETS
474 /*
475 TString arrname="fCandidateArray";
476 TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
477 fCandidateArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));
478 if (!fCandidateArray) {
479 Printf("%s array not found",candarrname.Data());
480 InputEvent()->GetList()->ls();
481 return;
482 }
483 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
484
485 TString arrSBname="fSideBandArray";
486 TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
487 fSideBandArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));
488 if (fCandidateType==1 && !fSideBandArray) {
489 Printf("%s array not found",sbarrname.Data());
490 InputEvent()->GetList()->ls();
491 return;
492 }
493 //Printf("nSBCand or Bkg found %d",fSideBandArray->GetEntriesFast());
494 */
495
496 //Histograms
497 TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
498 TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
499 TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
500 TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
501 TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
502 TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
503 TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
504 TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
505 TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
506 TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
507 TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
508 TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
509 TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
510 TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
511 TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
512 TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
513 THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
514
515 hstat->Fill(0);
516
517 // fix for temporary bug in ESDfilter
518 // the AODs with null vertex pointer didn't pass the PhysSel
519 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
520
521 //Event selection
522 //Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
523 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
524 if(!iseventselected) return;
525
526 hstat->Fill(1);
527
528 //retrieve charm candidates selected
529 Int_t candidates = fCandidateArray->GetEntriesFast();
530
531 //trigger on jets
532
533 Int_t njets=GetJetContainer()->GetNJets();
534 if(njets == 0) {
535 hstat->Fill(6, candidates);
536 hNDPerEvNoJet->Fill(candidates);
537 for(Int_t iD=0;iD<candidates;iD++){
538 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
539 if(!cand) continue;
540 hptDPerEvNoJet->Fill(cand->Pt());
541
542 }
543 return;
544
545 }
546
547 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
548
549 for(Int_t ic = 0; ic < candidates; ic++) {
550
551 // D* candidates
552 AliVParticle* charm=0x0;
553 charm=(AliVParticle*)fCandidateArray->At(ic);
554 if(!charm) continue;
555 hstat->Fill(2);
556
557 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
558
559 if(fCandidateType==kDstartoKpipi) {
560 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
561 Double_t deltamass= dstar->DeltaInvMass();
562 candsparse[3]=deltamass;
563 hnspDstandalone->Fill(candsparse);
564 }
565 if(fCandidateType==kD0toKpi){
566 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
567 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
568
569 Double_t masses[2];
570 //Int_t pdgdaughtersD0[2]={211,321};//pi,K
571 //Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
572
573 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
574 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
575 if(isselected==1 || isselected==3) {
576 candsparse[3]=masses[0];
577 hnspDstandalone->Fill(candsparse);
578 }
579 if(isselected>=2){
580 candsparse[3]=masses[1];
581 hnspDstandalone->Fill(candsparse);
582
583 }
584 }
585 }
586
587 // we start with jets
588 Double_t ejet = 0;
589 Double_t phiJet = 0;
590 Double_t etaJet = 0;
591 Double_t ptjet = 0;
592 Double_t leadingJet =0;
593
594 Int_t ntrarr=trackArr->GetEntriesFast();
595 hNtrArr->Fill(ntrarr);
596
597 for(Int_t i=0;i<ntrarr;i++){
598 AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
599 //reusing histograms
600 hPtJetTrks->Fill(jtrack->Pt());
601 hPhiJetTrks->Fill(jtrack->Phi());
602 hEtaJetTrks->Fill(jtrack->Eta());
603 hEjetTrks->Fill(jtrack->E());
604 }
605
606
607 //option to use only the leading jet
608 if(fLeadingJetOnly){
609 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
610 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
611 ptjet = jetL->Pt();
612 if(ptjet>leadingJet ) leadingJet = ptjet;
613
614 }
615 }
616
617 Int_t cntjet=0;
618 //loop over jets and consider only the leading jet in the event
619 for (Int_t iJets = 0; iJets<njets; iJets++) {
620 //Printf("Jet N %d",iJets);
621 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
622 //Printf("Corr task Accept Jet");
623 if(!AcceptJet(jet)) {
624 hstat->Fill(5);
625 continue;
626 }
627
628 //jets variables
629 ejet = jet->E();
630 phiJet = jet->Phi();
631 etaJet = jet->Eta();
632 ptjet = jet->Pt();
633
634 // choose the leading jet
635 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
636 //used for normalization
637 hstat->Fill(3);
638 cntjet++;
639 // fill energy, eta and phi of the jet
640 hEjet ->Fill(ejet);
641 hPhiJet ->Fill(phiJet);
642 hEtaJet ->Fill(etaJet);
643 hPtJet ->Fill(ptjet);
644
645 //loop on jet particles
646 Int_t ntrjet= jet->GetNumberOfTracks();
647 for(Int_t itrk=0;itrk<ntrjet;itrk++){
648
649 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
650 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
651
652 }//end loop on jet tracks
653
654 if(candidates==0){
655 hstat->Fill(7);
656 hPtJetPerEvNoD->Fill(jet->Pt());
657 }
658
659 //Printf("N candidates %d ", candidates);
660 for(Int_t ic = 0; ic < candidates; ic++) {
661
662 // D* candidates
663 AliVParticle* charm=0x0;
664 charm=(AliVParticle*)fCandidateArray->At(ic);
665 if(!charm) continue;
666
667
668 FlagFlavour(charm, jet);
669 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
670
671 FillHistogramsRecoJetCorr(charm, jet);
672
673 }
674 //retrieve side band background candidates for Dstar
675 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
676
677 for(Int_t ib=0;ib<nsbcand;ib++){
678 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
679 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
680 if(!sbcand) continue;
681 SideBandBackground(sbcand,jet);
682 }
683 if(fUseMCInfo){
684 AliAODRecoDecayHF* charmbg = 0x0;
685 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
686 if(!charmbg) continue;
687 MCBackground(charmbg,jet);
688 }
689 }
690 } // end of jet loop
691
692 hNJetPerEv->Fill(cntjet);
693 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
694 PostData(1,fmyOutput);
695 PostData(3,fmyOutputF);
696}
697
698//_______________________________________________________________________________
699
700void AliAnalysisTaskFlavourFilterAndJetCorrelations::Terminate(Option_t*)
701{
702 // The Terminate() function is the last function to be called during
703 // a query. It always runs on the client, it can be used to present
704 // the results graphically or save the results to file.
705
706 Info("Terminate"," terminate");
707 AliAnalysisTaskSE::Terminate();
708
709 fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
710 if (!fmyOutput) {
711 printf("ERROR: fmyOutput not available\n");
712 return;
713 }
714}
715
716//_______________________________________________________________________________
717
718void AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
719 Float_t mass=0;
720 Int_t abspdg=TMath::Abs(pdg);
721
722 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
723 // compute the Delta mass for the D*
724 if(fCandidateType==kDstartoKpipi){
725 Float_t mass1=0;
726 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
727 mass = mass-mass1;
728 }
729
730 fMinMass = mass-range;
731 fMaxMass = mass+range;
732
733 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
734 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
735}
736
737//_______________________________________________________________________________
738
739void AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
740 if(uplimit>lowlimit)
741 {
742 fMinMass = lowlimit;
743 fMaxMass = uplimit;
744 }
745 else{
746 printf("Error! Lower limit larger than upper limit!\n");
747 fMinMass = uplimit - uplimit*0.2;
748 fMaxMass = uplimit;
749 }
750}
751
752//_______________________________________________________________________________
753
754Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
755 if(nptbins>30) {
756 AliInfo("Maximum number of bins allowed is 30!");
757 return kFALSE;
758 }
759 if(!width) return kFALSE;
760 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
761 return kTRUE;
762}
763
764//_______________________________________________________________________________
765
766Double_t AliAnalysisTaskFlavourFilterAndJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
767 if(!part || !jet){
768 printf("Particle or jet do not exist!\n");
769 return -99;
770 }
771 Double_t p[3],pj[3];
772 Bool_t okpp=part->PxPyPz(p);
773 Bool_t okpj=jet->PxPyPz(pj);
774 if(!okpp || !okpj){
775 printf("Problems getting momenta\n");
776 return -999;
777 }
778 Double_t pjet=jet->P();
779 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
780 return z;
781}
782
783//_______________________________________________________________________________
784
785Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForAnalysis(){
786
787 // Statistics
788 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
789 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
790 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
791 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
792 hstat->GetXaxis()->SetBinLabel(4,"N jets");
793 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
794 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
795 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
796 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
797 hstat->SetNdivisions(1);
798 fmyOutput->Add(hstat);
799
800 const Int_t nbinsmass=200;
801 const Int_t nbinsptjet=500;
802 const Int_t nbinsptD=100;
803 const Int_t nbinsz=100;
804 const Int_t nbinsphi=200;
805 const Int_t nbinseta=100;
806
807 const Float_t ptjetlims[2]={0.,200.};
808 const Float_t ptDlims[2]={0.,50.};
809 const Float_t zlims[2]={0.,1.2};
810 const Float_t philims[2]={0.,6.3};
811 const Float_t etalims[2]={-1.5,1.5};
812
813 if(fCandidateType==kDstartoKpipi)
814 {
815
816 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
817 hDiffSideBand->SetStats(kTRUE);
818 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
819 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
820 hDiffSideBand->Sumw2();
821 fmyOutput->Add(hDiffSideBand);
822
823
824 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
825 hPtPion->SetStats(kTRUE);
826 hPtPion->GetXaxis()->SetTitle("GeV/c");
827 hPtPion->GetYaxis()->SetTitle("Entries");
828 hPtPion->Sumw2();
829 fmyOutput->Add(hPtPion);
830
831 }
832 // D related histograms
833 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
834 hNDPerEvNoJet->Sumw2();
835 fmyOutput->Add(hNDPerEvNoJet);
836
837 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
838 hptDPerEvNoJet->Sumw2();
839 fmyOutput->Add(hptDPerEvNoJet);
840
841 const Int_t nAxisD=4;
842 const Int_t nbinsSparseD[nAxisD]={nbinsphi,nbinseta,nbinsptD,nbinsmass};
843 const Double_t minSparseD[nAxisD] ={philims[0],etalims[0],ptDlims[0],fMinMass};
844 const Double_t maxSparseD[nAxisD] ={philims[1],etalims[1],ptDlims[1],fMaxMass};
845 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
846 hsDstandalone->Sumw2();
847
848 fmyOutput->Add(hsDstandalone);
849
850 // jet related fistograms
851
852 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
853 hEjetTrks->Sumw2();
854 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
855 hPhiJetTrks->Sumw2();
856 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
857 hEtaJetTrks->Sumw2();
858 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
859 hPtJetTrks->Sumw2();
860
861 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
862 hEjet->Sumw2();
863 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
864 hPhiJet->Sumw2();
865 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
866 hEtaJet->Sumw2();
867 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
868 hPtJet->Sumw2();
869
870 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]);
871 hPtJetWithD->Sumw2();
872 //for the MC this histogram is filled with the real background
873 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]);
874 hPtJetWithDsb->Sumw2();
875
876 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
877 hdeltaRJetTracks->Sumw2();
878
879 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
880 hNtrArr->Sumw2();
881
882 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
883 hNJetPerEv->Sumw2();
884
885 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);
886 hNJetPerEvNoD->Sumw2();
887
888 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
889 hPtJetPerEvNoD->Sumw2();
890
891 fmyOutput->Add(hEjetTrks);
892 fmyOutput->Add(hPhiJetTrks);
893 fmyOutput->Add(hEtaJetTrks);
894 fmyOutput->Add(hPtJetTrks);
895 fmyOutput->Add(hEjet);
896 fmyOutput->Add(hPhiJet);
897 fmyOutput->Add(hEtaJet);
898 fmyOutput->Add(hPtJet);
899 fmyOutput->Add(hPtJetWithD);
900 fmyOutput->Add(hPtJetWithDsb);
901 fmyOutput->Add(hdeltaRJetTracks);
902 fmyOutput->Add(hNtrArr);
903 fmyOutput->Add(hNJetPerEv);
904 fmyOutput->Add(hNJetPerEvNoD);
905 fmyOutput->Add(hPtJetPerEvNoD);
906
907 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
908 hDeltaRD->Sumw2();
909 fmyOutput->Add(hDeltaRD);
910
911 //background (side bands for the Dstar and like sign for D0)
912 fJetRadius=GetJetContainer(0)->GetJetRadius();
913 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]);
914 hInvMassptD->SetStats(kTRUE);
915 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
916 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
917 hInvMassptD->Sumw2();
918
919 fmyOutput->Add(hInvMassptD);
920
921 if(fUseMCInfo){
922 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]);
923 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
924 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
925 hInvMassptDbg->Sumw2();
926 fmyOutput->Add(hInvMassptDbg);
927
928 }
929 Double_t pi=TMath::Pi(), philims2[2];
930 philims2[0]=-pi*1./2.;
931 philims2[1]=pi*3./2.;
932 const Int_t nAxis=6;
933 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2};
934 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5};
935 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5};
936 THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, and mass", nAxis, nbinsSparse, minSparse, maxSparse);
937 hsDphiz->Sumw2();
938
939 fmyOutput->Add(hsDphiz);
940
941 PostData(1, fmyOutput);
942
943 return kTRUE;
944}
945
946//_______________________________________________________________________________
947
948Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForFiltering()
949{
950 //
951 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
952 //
953
954 // Statistics
955 TH1I* hstat = new TH1I("hstatF","Statistics",5,-0.5,4.5);
956 hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
957 hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
958 if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
959 else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
960 hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
961 if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
962 hstat->SetNdivisions(1);
963 fmyOutputF->Add(hstat);
964
965 TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
966 fmyOutputF->Add(hnCandEv);
967
968 // Invariant mass related histograms
969 const Int_t nbinsmass = 200;
970 TH2F *hInvMass = new TH2F("hInvMassptDF", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
971 hInvMass->SetStats(kTRUE);
972 hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
973 hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
974 fmyOutputF->Add(hInvMass);
975
976 if (fCandidateType==kDstartoKpipi) {
977 TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
978 fmyOutputF->Add(hnSBCandEv);
979
980 TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
981 hPtPion->SetStats(kTRUE);
982 hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
983 hPtPion->GetYaxis()->SetTitle("entries");
984 fmyOutputF->Add(hPtPion);
985 }
986 PostData(3, fmyOutputF);
987 return kTRUE;
988}
989
990//_______________________________________________________________________________
991
992void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
993
994 Double_t ptD=candidate->Pt();
995 Double_t ptjet=jet->Pt();
996 Double_t deltaR=DeltaR(candidate,jet);
997 Double_t deltaphi = jet->Phi()-candidate->Phi();
998 if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
999 if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
1000 Double_t z=Z(candidate,jet);
1001
1002 TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
1003 hDeltaRD->Fill(deltaR);
1004 if(fUseReco){
1005 if(fCandidateType==kD0toKpi) {
1006 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1007 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
1008
1009 }
1010
1011 if(fCandidateType==kDstartoKpipi) {
1012 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1013 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
1014
1015 }
1016 } else{ //generated level
1017 //AliInfo("Non reco");
1018 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
1019 }
1020}
1021
1022//_______________________________________________________________________________
1023
1024void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
1025
1026
1027 Double_t masses[2]={0.,0.};
1028 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1029 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1030
1031 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1032 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1033
1034 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1035 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1036 Double_t point[5]={z,dPhi,ptj,ptD,masses[0]};
1037
1038 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1039 if(isselected==1 || isselected==3) {
1040
1041 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
1042
1043 FillMassHistograms(masses[0], ptD, deltaR);
1044 hsDphiz->Fill(point,1.);
1045 }
1046 if(isselected>=2) {
1047 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
1048
1049 FillMassHistograms(masses[1], ptD, deltaR);
1050 point[4]=masses[1];
1051 hsDphiz->Fill(point,1.);
1052 }
1053
1054}
1055
1056//_______________________________________________________________________________
1057
1058void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
1059 //dPhi and z not used at the moment,but will be (re)added
1060
1061 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1062 Double_t deltamass= dstar->DeltaInvMass();
1063 //Double_t massD0= dstar->InvMassD0();
1064
1065
1066 TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
1067 hPtPion->Fill(softpi->Pt());
1068
1069 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1070 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1071 Int_t isSB=IsDzeroSideBand(dstar);
1072
1073 Double_t point[6]={z,dPhi,ptj,ptD,deltamass,isSB};
1074
1075 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
1076
1077 FillMassHistograms(deltamass, ptD, deltaR);
1078 hsDphiz->Fill(point,1.);
1079
1080}
1081
1082//_______________________________________________________________________________
1083
1084void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
1085
1086 Double_t pdgmass=0;
1087 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1088 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1089 Double_t point[6]={z,dPhi,ptjet,ptD,pdgmass,0};
1090
1091 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1092 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1093 point[4]=pdgmass;
1094
1095 if(deltaR<fJetRadius) {
1096 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1097 hsDphiz->Fill(point,1.);
1098 }
1099
1100}
1101
1102//_______________________________________________________________________________
1103
1104void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
1105
1106 if(deltaR<fJetRadius) {
1107 TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
1108 hInvMassptD->Fill(mass,ptD);
1109 }
1110}
1111
1112//________________________________________________________________________________
1113
1114void AliAnalysisTaskFlavourFilterAndJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
1115 Double_t deltaR=DeltaR(charm, jet);
1116 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1117 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1118 if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
1119
1120 return;
1121
1122}
1123
1124//_______________________________________________________________________________
1125
1126void AliAnalysisTaskFlavourFilterAndJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1127
1128 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1129 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1130 // has a width of ~5 sigmas. Two band needed for opening angle considerations
1131 TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
1132 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
1133
1134 Double_t deltaM=candDstar->DeltaInvMass();
1135 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1136 Double_t ptD=candDstar->Pt();
1137 Double_t ptjet=jet->Pt();
1138 Double_t dPhi=jet->Phi()-candDstar->Phi();
1139 Double_t deltaR=DeltaR(candDstar,jet);
1140 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1141 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1142
1143 Int_t isSideBand=IsDzeroSideBand(candDstar);
1144 //fill the background histogram with the side bands or when looking at MC with the real background
1145 if(isSideBand==1){
1146 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1147 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1148
1149 if(deltaR<fJetRadius){ // evaluate in the near side
1150 //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
1151 hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
1152 }
1153 }
1154}
1155
1156//_______________________________________________________________________________
1157
1158void AliAnalysisTaskFlavourFilterAndJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
1159
1160 //need mass, deltaR, pt jet, ptD
1161 //two cases: D0 and Dstar
1162
1163 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1164 if(!isselected) return;
1165
1166 Double_t ptD=candbg->Pt();
1167 Double_t ptjet=jet->Pt();
1168 Double_t deltaR=DeltaR(candbg,jet);
1169
1170 TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
1171 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
1172
1173
1174 if(fCandidateType==kDstartoKpipi){
1175 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1176 Double_t deltaM=dstarbg->DeltaInvMass();
1177 hInvMassptDbg->Fill(deltaM,ptD);
1178 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
1179 }
1180
1181 if(fCandidateType==kD0toKpi){
1182 Double_t masses[2]={0.,0.};
1183 Int_t pdgdaughtersD0[2]={211,321};//pi,K
1184 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1185
1186 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1187 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1188
1189
1190 if(isselected==1 || isselected==3) {
1191 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
1192 hInvMassptDbg->Fill(masses[0],ptD);
1193 }
1194 if(isselected>=2) {
1195 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
1196 hInvMassptDbg->Fill(masses[1],ptD);
1197 }
1198
1199
1200 }
1201}
1202
1203//_______________________________________________________________________________
1204
1205Float_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1206 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1207
1208 if(!p1 || !p2) return -1;
1209 Double_t phi1=p1->Phi(),eta1=p1->Eta();
1210 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1211
1212 Double_t dPhi=phi1-phi2;
1213 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1214 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1215
1216 Double_t dEta=eta1-eta2;
1217 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1218 return deltaR;
1219
1220}
1221
1222//_______________________________________________________________________________
1223
1224Int_t AliAnalysisTaskFlavourFilterAndJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1225
1226 Double_t ptD=candDstar->Pt();
1227 Int_t bin = fCuts->PtBin(ptD);
1228 if (bin < 0){
1229 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1230 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1231 return -1; // out of bounds
1232 }
1233
1234 Double_t invM=candDstar->InvMassD0();
1235 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1236
1237 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1238 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1239
1240 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1241 else return 0;
1242
1243}