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