]> 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
implementing the VZERO selection in p
[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>
33
34#include "AliAnalysisTaskFlavourJetCorrelations.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(AliAnalysisTaskFlavourJetCorrelations)
51
52
53//_______________________________________________________________________________
54
55AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
56AliAnalysisTaskEmcalJet("",kFALSE),
57fUseMCInfo(kTRUE),
58fUseReco(kTRUE),
59fCandidateType(),
60fPDGmother(),
61fNProngs(),
62fPDGdaughters(),
63fBranchName(),
931f8b00 64fCuts(0),
65fMinMass(),
66fMaxMass(),
67fJetArrName(0),
68fCandArrName(0),
69fLeadingJetOnly(kFALSE),
bbb94467 70fJetRadius(0),
71fCandidateArray(0),
72fSideBandArray(0)
931f8b00 73{
74 //
75 // Default ctor
76 //
77 //SetMakeGeneralHistograms(kTRUE);
78
79}
80
81//_______________________________________________________________________________
82
83AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
84AliAnalysisTaskEmcalJet(name,kFALSE),
85fUseMCInfo(kTRUE),
86fUseReco(kTRUE),
87fCandidateType(),
88fPDGmother(),
89fNProngs(),
90fPDGdaughters(),
91fBranchName(),
931f8b00 92fCuts(0),
93fMinMass(),
94fMaxMass(),
95fJetArrName(0),
96fCandArrName(0),
97fLeadingJetOnly(kFALSE),
bbb94467 98fJetRadius(0),
99fCandidateArray(0),
100fSideBandArray(0)
931f8b00 101{
102 //
103 // Constructor. Initialization of Inputs and Outputs
104 //
105
106 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
107 fCuts=cuts;
108 fCandidateType=candtype;
109 const Int_t nptbins=fCuts->GetNPtBins();
110 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};
111
112 switch(fCandidateType){
113 case 0:
114 fPDGmother=421;
115 fNProngs=2;
116 fPDGdaughters[0]=211;//pi
117 fPDGdaughters[1]=321;//K
118 fPDGdaughters[2]=0; //empty
119 fPDGdaughters[3]=0; //empty
120 fBranchName="D0toKpi";
121 fCandArrName="D0";
122 break;
123 case 1:
124 fPDGmother=413;
125 fNProngs=3;
126 fPDGdaughters[1]=211;//pi soft
127 fPDGdaughters[0]=421;//D0
128 fPDGdaughters[2]=211;//pi fromD0
129 fPDGdaughters[3]=321; // K from D0
130 fBranchName="Dstar";
131 fCandArrName="Dstar";
132
133 if(nptbins<=13){
134 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
135 } else {
136 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
137 }
138 break;
139 default:
140 printf("%d not accepted!!\n",fCandidateType);
141 break;
142 }
143
144 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
145 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
146
bbb94467 147 DefineInput(1, TClonesArray::Class());
148 DefineInput(2, TClonesArray::Class());
149
931f8b00 150 DefineOutput(1,TList::Class()); // histos
151 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
152
153}
154
155//_______________________________________________________________________________
156
157AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
158 //
159 // destructor
160 //
161
162 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
163
931f8b00 164 delete fCuts;
165
166}
167
168//_______________________________________________________________________________
169
170void AliAnalysisTaskFlavourJetCorrelations::Init(){
171 //
172 // Initialization
173 //
174
175 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
176 switch(fCandidateType){
177 case 0:
178 {
179 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
180 copyfCuts->SetName("AnalysisCutsDzero");
181 // Post the data
182 PostData(2,copyfCuts);
183 }
184 break;
185 case 1:
186 {
187 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
188 copyfCuts->SetName("AnalysisCutsDStar");
189 // Post the cuts
190 PostData(2,copyfCuts);
191 }
192 break;
193 default:
194 return;
195 }
196
197 return;
198}
199
200//_______________________________________________________________________________
201
202void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
203 // output
204 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
76bf81f2 205 AliAnalysisTaskEmcal::UserCreateOutputObjects();
931f8b00 206 // define histograms
76bf81f2 207 // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
931f8b00 208 DefineHistoForAnalysis();
76bf81f2 209 PostData(1,fOutput);
931f8b00 210
211 return;
212}
213
214//_______________________________________________________________________________
215
b5d0ba13 216Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
931f8b00 217{
76bf81f2 218 // user exec from AliAnalysisTaskEmcal is used
219
931f8b00 220 // Load the event
221 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
222
223 TClonesArray *arrayDStartoD0pi=0;
224 TClonesArray *trackArr = 0;
931f8b00 225
226 if (!aodEvent && AODEvent() && IsStandardAOD()) {
227
228 // In case there is an AOD handler writing a standard AOD, use the AOD
229 // event in memory rather than the input (ESD) event.
230 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
231
232 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
233 // have to taken from the AOD event hold by the AliAODExtension
234 AliAODHandler* aodHandler = (AliAODHandler*)
235 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
236 if(aodHandler->GetExtensions()) {
237 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
238 AliAODEvent *aodFromExt = ext->GetAOD();
239 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
240 }
241 } else {
242 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
243 }
244
245 if (!arrayDStartoD0pi) {
246 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
247 // return;
248 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
249
250 TClonesArray* mcArray = 0x0;
251 if (fUseMCInfo) {
252 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
253 if (!mcArray) {
254 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
b5d0ba13 255 return kFALSE;
931f8b00 256 }
257 }
258
259 //retrieve jets
260 trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
261 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
262 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
263 fJetRadius=GetJetContainer(0)->GetJetRadius();
264
265 if(!trackArr){
266 AliInfo(Form("Could not find the track array\n"));
b5d0ba13 267 return kFALSE;
931f8b00 268 }
269
bbb94467 270
c3c1ab31 271 fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
b5d0ba13 272 if (!fCandidateArray) return kFALSE;
c3c1ab31 273 if (fCandidateType==1) {
274 fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
b5d0ba13 275 if (!fSideBandArray) return kFALSE;
c3c1ab31 276 }
277 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
931f8b00 278
279 //Histograms
76bf81f2 280 TH1I* hstat=(TH1I*)fOutput->FindObject("hstat");
281 TH1F* hPtJetTrks=(TH1F*)fOutput->FindObject("hPtJetTrks");
282 TH1F* hPhiJetTrks=(TH1F*)fOutput->FindObject("hPhiJetTrks");
283 TH1F* hEtaJetTrks=(TH1F*)fOutput->FindObject("hEtaJetTrks");
284 TH1F* hEjetTrks=(TH1F*)fOutput->FindObject("hEjetTrks");
285 TH1F* hPtJet=(TH1F*)fOutput->FindObject("hPtJet");
286 TH1F* hPhiJet=(TH1F*)fOutput->FindObject("hPhiJet");
287 TH1F* hEtaJet=(TH1F*)fOutput->FindObject("hEtaJet");
288 TH1F* hEjet=(TH1F*)fOutput->FindObject("hEjet");
289 TH1F* hNtrArr=(TH1F*)fOutput->FindObject("hNtrArr");
290 TH1F* hNJetPerEv=(TH1F*)fOutput->FindObject("hNJetPerEv");
291 TH1F* hdeltaRJetTracks=(TH1F*)fOutput->FindObject("hdeltaRJetTracks");
292 TH1F* hNDPerEvNoJet=(TH1F*)fOutput->FindObject("hNDPerEvNoJet");
293 TH1F* hptDPerEvNoJet=(TH1F*)fOutput->FindObject("hptDPerEvNoJet");
294 TH1F* hNJetPerEvNoD=(TH1F*)fOutput->FindObject("hNJetPerEvNoD");
295 TH1F* hPtJetPerEvNoD=(TH1F*)fOutput->FindObject("hPtJetPerEvNoD");
296 THnSparseF* hnspDstandalone=(THnSparseF*)fOutput->FindObject("hsDstandalone");
d9126a4c 297
931f8b00 298 hstat->Fill(0);
299
300 // fix for temporary bug in ESDfilter
301 // the AODs with null vertex pointer didn't pass the PhysSel
b5d0ba13 302 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
931f8b00 303
304 //Event selection
305 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
306 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
b5d0ba13 307 if(!iseventselected) return kFALSE;
931f8b00 308
309 hstat->Fill(1);
931f8b00 310
311 //retrieve charm candidates selected
bbb94467 312 Int_t candidates = fCandidateArray->GetEntriesFast();
d9126a4c
CB		 313
314 //trigger on jets
931f8b00 315
d9126a4c
CB		 316 Int_t njets=GetJetContainer()->GetNJets();
317 if(njets == 0) {
318 hstat->Fill(6, candidates);
319 hNDPerEvNoJet->Fill(candidates);
320 for(Int_t iD=0;iD<candidates;iD++){
bbb94467 321 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
d9126a4c
CB		 322 if(!cand) continue;
323 hptDPerEvNoJet->Fill(cand->Pt());
324
325 }
b5d0ba13 326 return kFALSE;
d9126a4c
CB		 327
328 }
1353fa73 329
330 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
331
332 for(Int_t ic = 0; ic < candidates; ic++) {
333
334 // D* candidates
bbb94467 335 AliAODRecoDecayHF* charm=0x0;
336 AliAODRecoCascadeHF* dstar=0x0;
337
338
339 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
1353fa73 340 if(!charm) continue;
bbb94467 341
342 if(fCandidateType==kDstartoKpipi){
343 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
344 if(!dstar) continue;
345 }
346
1353fa73 347 hstat->Fill(2);
348
349 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
350
351 if(fCandidateType==kDstartoKpipi) {
8577f7bc 352 if(fUseReco){
353 Double_t deltamass= dstar->DeltaInvMass();
354 candsparse[3]=deltamass;
355 } else candsparse[3] = 0.145; //for generated
1353fa73 356 hnspDstandalone->Fill(candsparse);
357 }
358 if(fCandidateType==kD0toKpi){
8577f7bc 359 if(fUseReco){
360 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
361 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
362 //not a further selection,just to get the value of isselected for the D0
363 Double_t masses[2];
364 Int_t pdgdaughtersD0[2]={211,321};//pi,K
365 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
366
367 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
368 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
369 if(isselected==1 || isselected==3) {
370 candsparse[3]=masses[0];
371 hnspDstandalone->Fill(candsparse);
372 }
373 if(isselected>=2){
374 candsparse[3]=masses[1];
375 hnspDstandalone->Fill(candsparse);
376
377 }
378 } else { //generated
379 Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
380 candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1353fa73 381 hnspDstandalone->Fill(candsparse);
1353fa73 382 }
383 }
384 }
385
931f8b00 386 // we start with jets
387 Double_t ejet = 0;
388 Double_t phiJet = 0;
389 Double_t etaJet = 0;
390 Double_t ptjet = 0;
391 Double_t leadingJet =0;
392
393 Int_t ntrarr=trackArr->GetEntriesFast();
394 hNtrArr->Fill(ntrarr);
395
396 for(Int_t i=0;i<ntrarr;i++){
397 AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
398 //reusing histograms
399 hPtJetTrks->Fill(jtrack->Pt());
400 hPhiJetTrks->Fill(jtrack->Phi());
401 hEtaJetTrks->Fill(jtrack->Eta());
402 hEjetTrks->Fill(jtrack->E());
403 }
404
405
406 //option to use only the leading jet
407 if(fLeadingJetOnly){
408 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
409 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
410 ptjet = jetL->Pt();
411 if(ptjet>leadingJet ) leadingJet = ptjet;
412
413 }
414 }
415
416 Int_t cntjet=0;
417 //loop over jets and consider only the leading jet in the event
418 for (Int_t iJets = 0; iJets<njets; iJets++) {
419 //Printf("Jet N %d",iJets);
420 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
421 //Printf("Corr task Accept Jet");
422 if(!AcceptJet(jet)) {
423 hstat->Fill(5);
424 continue;
425 }
426
427 //jets variables
428 ejet = jet->E();
429 phiJet = jet->Phi();
430 etaJet = jet->Eta();
431 ptjet = jet->Pt();
432
433 // choose the leading jet
434 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
435 //used for normalization
436 hstat->Fill(3);
437 cntjet++;
438 // fill energy, eta and phi of the jet
439 hEjet ->Fill(ejet);
440 hPhiJet ->Fill(phiJet);
441 hEtaJet ->Fill(etaJet);
442 hPtJet ->Fill(ptjet);
443
444 //loop on jet particles
445 Int_t ntrjet= jet->GetNumberOfTracks();
446 for(Int_t itrk=0;itrk<ntrjet;itrk++){
447
448 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
449 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
450
451 }//end loop on jet tracks
452
d9126a4c
CB		 453 if(candidates==0){
454 hstat->Fill(7);
455 hPtJetPerEvNoD->Fill(jet->Pt());
456 }
457
931f8b00 458 //Printf("N candidates %d ", candidates);
459 for(Int_t ic = 0; ic < candidates; ic++) {
460
461 // D* candidates
462 AliVParticle* charm=0x0;
bbb94467 463 charm=(AliVParticle*)fCandidateArray->At(ic);
931f8b00 464 if(!charm) continue;
1353fa73 465
931f8b00 466
467 FlagFlavour(charm, jet);
468 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
469
8577f7bc 470 FillHistogramsRecoJetCorr(charm, jet, aodEvent);
931f8b00 471
472 }
473 //retrieve side band background candidates for Dstar
bbb94467 474 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
931f8b00 475
476 for(Int_t ib=0;ib<nsbcand;ib++){
477 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
bbb94467 478 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
931f8b00 479 if(!sbcand) continue;
480 SideBandBackground(sbcand,jet);
481 }
482 if(fUseMCInfo){
483 AliAODRecoDecayHF* charmbg = 0x0;
bbb94467 484 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
931f8b00 485 if(!charmbg) continue;
486 MCBackground(charmbg,jet);
487 }
488 }
489 } // end of jet loop
490
491 hNJetPerEv->Fill(cntjet);
d9126a4c 492 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
76bf81f2 493 PostData(1,fOutput);
b5d0ba13 494 return kTRUE;
931f8b00 495}
496
497//_______________________________________________________________________________
498
499void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
500{
501 // The Terminate() function is the last function to be called during
502 // a query. It always runs on the client, it can be used to present
503 // the results graphically or save the results to file.
504
505 Info("Terminate"," terminate");
506 AliAnalysisTaskSE::Terminate();
507
76bf81f2 508 fOutput = dynamic_cast<TList*> (GetOutputData(1));
509 if (!fOutput) {
510 printf("ERROR: fOutput not available\n");
931f8b00 511 return;
512 }
513}
514
515//_______________________________________________________________________________
516
517void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
518 Float_t mass=0;
519 Int_t abspdg=TMath::Abs(pdg);
520
521 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
522 // compute the Delta mass for the D*
523 if(fCandidateType==kDstartoKpipi){
524 Float_t mass1=0;
525 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
526 mass = mass-mass1;
527 }
528
529 fMinMass = mass-range;
530 fMaxMass = mass+range;
531
532 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
533 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
534}
535
536//_______________________________________________________________________________
537
538void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
539 if(uplimit>lowlimit)
540 {
541 fMinMass = lowlimit;
542 fMaxMass = uplimit;
543 }
544 else{
545 printf("Error! Lower limit larger than upper limit!\n");
546 fMinMass = uplimit - uplimit*0.2;
547 fMaxMass = uplimit;
548 }
549}
550
551//_______________________________________________________________________________
552
553Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
554 if(nptbins>30) {
555 AliInfo("Maximum number of bins allowed is 30!");
556 return kFALSE;
557 }
558 if(!width) return kFALSE;
559 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
560 return kTRUE;
561}
562
563//_______________________________________________________________________________
564
565Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
566 if(!part || !jet){
567 printf("Particle or jet do not exist!\n");
568 return -99;
569 }
570 Double_t p[3],pj[3];
571 Bool_t okpp=part->PxPyPz(p);
572 Bool_t okpj=jet->PxPyPz(pj);
573 if(!okpp || !okpj){
574 printf("Problems getting momenta\n");
575 return -999;
576 }
577 Double_t pjet=jet->P();
578 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
579 return z;
580}
581
582//_______________________________________________________________________________
583
584Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
585
586 // Statistics
d9126a4c 587 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
931f8b00 588 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
589 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
d9126a4c 590 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
931f8b00 591 hstat->GetXaxis()->SetBinLabel(4,"N jets");
592 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
593 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
d9126a4c
CB		 594 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
595 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
931f8b00 596 hstat->SetNdivisions(1);
76bf81f2 597 fOutput->Add(hstat);
931f8b00 598
599 const Int_t nbinsmass=200;
600 const Int_t nbinsptjet=500;
601 const Int_t nbinsptD=100;
602 const Int_t nbinsz=100;
603 const Int_t nbinsphi=200;
1353fa73 604 const Int_t nbinseta=100;
931f8b00 605
606 const Float_t ptjetlims[2]={0.,200.};
607 const Float_t ptDlims[2]={0.,50.};
608 const Float_t zlims[2]={0.,1.2};
609 const Float_t philims[2]={0.,6.3};
1353fa73 610 const Float_t etalims[2]={-1.5,1.5};
931f8b00 611
612 if(fCandidateType==kDstartoKpipi)
613 {
614
615 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
616 hDiffSideBand->SetStats(kTRUE);
617 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
618 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
619 hDiffSideBand->Sumw2();
76bf81f2 620 fOutput->Add(hDiffSideBand);
931f8b00 621
622
623 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
624 hPtPion->SetStats(kTRUE);
625 hPtPion->GetXaxis()->SetTitle("GeV/c");
626 hPtPion->GetYaxis()->SetTitle("Entries");
627 hPtPion->Sumw2();
76bf81f2 628 fOutput->Add(hPtPion);
931f8b00 629
630 }
d9126a4c 631 // D related histograms
1353fa73 632 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
633 hNDPerEvNoJet->Sumw2();
76bf81f2 634 fOutput->Add(hNDPerEvNoJet);
1353fa73 635
636 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
637 hptDPerEvNoJet->Sumw2();
76bf81f2 638 fOutput->Add(hptDPerEvNoJet);
1353fa73 639
640 const Int_t nAxisD=4;
e3017b3f 641 const Int_t nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
642 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
643 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
1353fa73 644 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
645 hsDstandalone->Sumw2();
646
76bf81f2 647 fOutput->Add(hsDstandalone);
1353fa73 648
931f8b00 649 // jet related fistograms
650
651 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
652 hEjetTrks->Sumw2();
653 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
654 hPhiJetTrks->Sumw2();
1353fa73 655 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 656 hEtaJetTrks->Sumw2();
657 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
658 hPtJetTrks->Sumw2();
659
660 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
661 hEjet->Sumw2();
662 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
663 hPhiJet->Sumw2();
1353fa73 664 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 665 hEtaJet->Sumw2();
666 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
667 hPtJet->Sumw2();
668
669 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]);
670 hPtJetWithD->Sumw2();
671 //for the MC this histogram is filled with the real background
672 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]);
673 hPtJetWithDsb->Sumw2();
674
675 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
676 hdeltaRJetTracks->Sumw2();
677
678 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
679 hNtrArr->Sumw2();
d9126a4c 680
931f8b00 681 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
682 hNJetPerEv->Sumw2();
683
d9126a4c
CB		 684 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);
685 hNJetPerEvNoD->Sumw2();
686
687 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
688 hPtJetPerEvNoD->Sumw2();
689
76bf81f2 690 fOutput->Add(hEjetTrks);
691 fOutput->Add(hPhiJetTrks);
692 fOutput->Add(hEtaJetTrks);
693 fOutput->Add(hPtJetTrks);
694 fOutput->Add(hEjet);
695 fOutput->Add(hPhiJet);
696 fOutput->Add(hEtaJet);
697 fOutput->Add(hPtJet);
698 fOutput->Add(hPtJetWithD);
699 fOutput->Add(hPtJetWithDsb);
700 fOutput->Add(hdeltaRJetTracks);
701 fOutput->Add(hNtrArr);
702 fOutput->Add(hNJetPerEv);
703 fOutput->Add(hNJetPerEvNoD);
704 fOutput->Add(hPtJetPerEvNoD);
931f8b00 705
706 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
707 hDeltaRD->Sumw2();
76bf81f2 708 fOutput->Add(hDeltaRD);
931f8b00 709
710 //background (side bands for the Dstar and like sign for D0)
711 fJetRadius=GetJetContainer(0)->GetJetRadius();
712 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]);
713 hInvMassptD->SetStats(kTRUE);
714 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
715 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
716 hInvMassptD->Sumw2();
717
76bf81f2 718 fOutput->Add(hInvMassptD);
931f8b00 719
720 if(fUseMCInfo){
721 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]);
722 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
723 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
724 hInvMassptDbg->Sumw2();
76bf81f2 725 fOutput->Add(hInvMassptDbg);
931f8b00 726
727 }
728 Double_t pi=TMath::Pi(), philims2[2];
729 philims2[0]=-pi*1./2.;
730 philims2[1]=pi*3./2.;
ed5a88ff 731 const Int_t nAxis=7;
732 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
733 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
734 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
735 THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. For SB or not and within the jet cone or not", nAxis, nbinsSparse, minSparse, maxSparse);
931f8b00 736 hsDphiz->Sumw2();
737
76bf81f2 738 fOutput->Add(hsDphiz);
931f8b00 739
76bf81f2 740 PostData(1, fOutput);
931f8b00 741
742 return kTRUE;
743}
744
745//_______________________________________________________________________________
746
8577f7bc 747void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet, AliAODEvent* aodEvent){
931f8b00 748
749 Double_t ptD=candidate->Pt();
750 Double_t ptjet=jet->Pt();
751 Double_t deltaR=DeltaR(candidate,jet);
752 Double_t deltaphi = jet->Phi()-candidate->Phi();
753 if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
754 if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
755 Double_t z=Z(candidate,jet);
756
76bf81f2 757 TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
931f8b00 758 hDeltaRD->Fill(deltaR);
759 if(fUseReco){
760 if(fCandidateType==kD0toKpi) {
761 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
8577f7bc 762
763 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, aodEvent);
931f8b00 764
765 }
766
767 if(fCandidateType==kDstartoKpipi) {
768 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
769 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
770
771 }
772 } else{ //generated level
773 //AliInfo("Non reco");
774 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
775 }
776}
777
778//_______________________________________________________________________________
779
780void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
781
931f8b00 782
783 Double_t masses[2]={0.,0.};
784 Int_t pdgdaughtersD0[2]={211,321};//pi,K
785 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
786
787 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
788 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
789
76bf81f2 790 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
791 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
ed5a88ff 792 Double_t point[7]={z,dPhi,ptj,ptD,masses[0],0, deltaR<fJetRadius ? 1 : 0};
8577f7bc 793 Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
931f8b00 794 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
795 if(isselected==1 || isselected==3) {
796
797 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
798
799 FillMassHistograms(masses[0], ptD, deltaR);
800 hsDphiz->Fill(point,1.);
801 }
802 if(isselected>=2) {
803 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
804
805 FillMassHistograms(masses[1], ptD, deltaR);
806 point[4]=masses[1];
807 hsDphiz->Fill(point,1.);
808 }
809
810}
811
812//_______________________________________________________________________________
813
814void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
815 //dPhi and z not used at the moment,but will be (re)added
816
817 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
818 Double_t deltamass= dstar->DeltaInvMass();
819 //Double_t massD0= dstar->InvMassD0();
820
821
76bf81f2 822 TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
931f8b00 823 hPtPion->Fill(softpi->Pt());
824
76bf81f2 825 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
826 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
ed5a88ff 827 Int_t isSB=0;//IsDzeroSideBand(dstar);
931f8b00 828
ed5a88ff 829 Double_t point[7]={z,dPhi,ptj,ptD,deltamass,isSB, deltaR<fJetRadius ? 1 : 0};
931f8b00 830
831 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
832
833 FillMassHistograms(deltamass, ptD, deltaR);
834 hsDphiz->Fill(point,1.);
835
836}
837
838//_______________________________________________________________________________
839
840void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
841
842 Double_t pdgmass=0;
76bf81f2 843 TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
844 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
ed5a88ff 845 Double_t point[7]={z,dPhi,ptjet,ptD,pdgmass,0, deltaR<fJetRadius ? 1 : 0};
931f8b00 846
847 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
848 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
849 point[4]=pdgmass;
8577f7bc 850 hsDphiz->Fill(point,1.);
931f8b00 851 if(deltaR<fJetRadius) {
852 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
931f8b00 853 }
854
855}
856
857//_______________________________________________________________________________
858
859void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
860
861 if(deltaR<fJetRadius) {
76bf81f2 862 TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
931f8b00 863 hInvMassptD->Fill(mass,ptD);
864 }
865}
866
1353fa73 867//________________________________________________________________________________
868
931f8b00 869void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
870 Double_t deltaR=DeltaR(charm, jet);
871 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
872 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
873 if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
874
875 return;
876
877}
878
879//_______________________________________________________________________________
880
881void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
882
883 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
884 // (expected detector resolution) on the left and right frm the D0 mass. Each band
885 // has a width of ~5 sigmas. Two band needed for opening angle considerations
76bf81f2 886 TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
887 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
888 THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
931f8b00 889
890 Double_t deltaM=candDstar->DeltaInvMass();
891 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
ed5a88ff 892 Double_t z=Z(candDstar,jet);
931f8b00 893 Double_t ptD=candDstar->Pt();
894 Double_t ptjet=jet->Pt();
895 Double_t dPhi=jet->Phi()-candDstar->Phi();
896 Double_t deltaR=DeltaR(candDstar,jet);
897 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
898 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
899
ed5a88ff 900 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
901 Double_t point[7]={z,dPhi,ptjet,ptD,deltaM,isSideBand, deltaR<fJetRadius ? 1 : 0};
931f8b00 902 //fill the background histogram with the side bands or when looking at MC with the real background
903 if(isSideBand==1){
904 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
905 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
ed5a88ff 906 hsDphiz->Fill(point,1.);
931f8b00 907 if(deltaR<fJetRadius){ // evaluate in the near side
908 //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
909 hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
910 }
911 }
912}
913
914//_______________________________________________________________________________
915
916void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
917
918 //need mass, deltaR, pt jet, ptD
919 //two cases: D0 and Dstar
920
921 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
922 if(!isselected) return;
923
924 Double_t ptD=candbg->Pt();
925 Double_t ptjet=jet->Pt();
926 Double_t deltaR=DeltaR(candbg,jet);
927
76bf81f2 928 TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
929 TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
931f8b00 930
931
932 if(fCandidateType==kDstartoKpipi){
933 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
934 Double_t deltaM=dstarbg->DeltaInvMass();
935 hInvMassptDbg->Fill(deltaM,ptD);
936 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
937 }
938
939 if(fCandidateType==kD0toKpi){
940 Double_t masses[2]={0.,0.};
941 Int_t pdgdaughtersD0[2]={211,321};//pi,K
942 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
943
944 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
945 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
946
947
948 if(isselected==1 || isselected==3) {
949 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
950 hInvMassptDbg->Fill(masses[0],ptD);
951 }
952 if(isselected>=2) {
953 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
954 hInvMassptDbg->Fill(masses[1],ptD);
955 }
956
957
958 }
959}
960
961//_______________________________________________________________________________
962
963Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
964 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
965
966 if(!p1 || !p2) return -1;
967 Double_t phi1=p1->Phi(),eta1=p1->Eta();
968 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
969
970 Double_t dPhi=phi1-phi2;
971 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
972 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
973
974 Double_t dEta=eta1-eta2;
975 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
976 return deltaR;
977
978}
979
980//_______________________________________________________________________________
981
982Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
983
984 Double_t ptD=candDstar->Pt();
985 Int_t bin = fCuts->PtBin(ptD);
986 if (bin < 0){
987 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
988 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
989 return -1; // out of bounds
990 }
991
992 Double_t invM=candDstar->InvMassD0();
993 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
994
995 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
996 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
997
998 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
999 else return 0;
1000
1001}