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