]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
used computed res if available.
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSED0Correlations.cxx
CommitLineData
a2ad7da1 1/**************************************************************************
2 * Copyright(c) 1998-2012, 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/* $Id$ */
17
18/////////////////////////////////////////////////////////////
19//
20// AliAnalysisTaskSE for D0 candidates (2Prongs)
21// and hadrons correlations
22//
23// Authors: Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
24// Fabio Colamaria, fabio.colamaria@ba.infn.it (correlations)
25/////////////////////////////////////////////////////////////
26
27#include <Riostream.h>
28#include <TClonesArray.h>
29#include <TCanvas.h>
30#include <TNtuple.h>
31#include <TTree.h>
32#include <TList.h>
33#include <TH1F.h>
34#include <TH2F.h>
35#include <THnSparse.h>
36#include <TDatabasePDG.h>
37
38#include <AliAnalysisDataSlot.h>
39#include <AliAnalysisDataContainer.h>
40#include "AliAnalysisManager.h"
41#include "AliESDtrack.h"
42#include "AliVertexerTracks.h"
43#include "AliAODHandler.h"
44#include "AliAODEvent.h"
45#include "AliAODVertex.h"
46#include "AliAODTrack.h"
47#include "AliAODMCHeader.h"
48#include "AliAODMCParticle.h"
49#include "AliAODRecoDecayHF2Prong.h"
50#include "AliAODRecoCascadeHF.h"
51#include "AliAnalysisVertexingHF.h"
52#include "AliAnalysisTaskSE.h"
53#include "AliAnalysisTaskSED0Correlations.h"
54#include "AliNormalizationCounter.h"
55
56using std::cout;
57using std::endl;
58
59ClassImp(AliAnalysisTaskSED0Correlations)
60
61
62//________________________________________________________________________
63AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations():
64AliAnalysisTaskSE(),
65 fNPtBinsCorr(0),
66 fBinLimsCorr(),
67 fPtThreshLow(),
68 fPtThreshUp(),
69 fEvents(0),
70 fAlreadyFilled(kFALSE),
71 fOutputMass(0),
72 fOutputCorr(0),
73 fOutputStudy(0),
74 fNentries(0),
75 fCutsD0(0),
76 fCutsTracks(0),
77 fReadMC(0),
78 fCounter(0),
79 fNPtBins(1),
80 fFillOnlyD0D0bar(0),
81 fIsSelectedCandidate(0),
82 fSys(0),
83 fIsRejectSDDClusters(0),
84 fFillGlobal(kTRUE)
85{
86 // Default constructor
87
88}
89
90//________________________________________________________________________
91AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *name,AliRDHFCutsD0toKpi* cutsD0, AliHFAssociatedTrackCuts* cutsTrk):
92 AliAnalysisTaskSE(name),
93 fNPtBinsCorr(0),
94 fBinLimsCorr(),
95 fPtThreshLow(),
96 fPtThreshUp(),
97 fEvents(0),
98 fAlreadyFilled(kFALSE),
99 fOutputMass(0),
100 fOutputCorr(0),
101 fOutputStudy(0),
102 fNentries(0),
103 fCutsD0(0),
104 fCutsTracks(cutsTrk),
105 fReadMC(0),
106 fCounter(0),
107 fNPtBins(1),
108 fFillOnlyD0D0bar(0),
109 fIsSelectedCandidate(0),
110 fSys(0),
111 fIsRejectSDDClusters(0),
112 fFillGlobal(kTRUE)
113{
114 // Default constructor
115
116 fNPtBins=cutsD0->GetNPtBins();
117
118 fCutsD0=cutsD0;
119
120 // Output slot #1 writes into a TList container (mass with cuts)
121 DefineOutput(1,TList::Class()); //My private output
122 // Output slot #2 writes into a TH1F container (number of events)
123 DefineOutput(2,TH1F::Class()); //My private output
124 // Output slot #3 writes into a AliRDHFD0toKpi container (cuts)
125 DefineOutput(3,AliRDHFCutsD0toKpi::Class()); //My private output
126 // Output slot #4 writes Normalization Counter
127 DefineOutput(4,AliNormalizationCounter::Class());
128 // Output slot #5 writes into a TList container (correl output)
129 DefineOutput(5,TList::Class()); //My private output
130 // Output slot #6 writes into a TList container (correl advanced)
131 DefineOutput(6,TList::Class()); //My private output
132 // Output slot #7 writes into a AliHFAssociatedTrackCuts container (cuts)
133 DefineOutput(7,AliHFAssociatedTrackCuts::Class()); //My private output
134}
135
136//________________________________________________________________________
137AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalysisTaskSED0Correlations &source):
138 AliAnalysisTaskSE(source),
139 fNPtBinsCorr(source.fNPtBinsCorr),
140 fBinLimsCorr(source.fBinLimsCorr),
141 fPtThreshLow(source.fPtThreshLow),
142 fPtThreshUp(source.fPtThreshUp),
143 fEvents(source.fEvents),
144 fAlreadyFilled(source.fAlreadyFilled),
145 fOutputMass(source.fOutputMass),
146 fOutputCorr(source.fOutputCorr),
147 fOutputStudy(source.fOutputStudy),
148 fNentries(source.fNentries),
149 fCutsD0(source.fCutsD0),
150 fCutsTracks(source.fCutsTracks),
151 fReadMC(source.fReadMC),
152 fCounter(source.fCounter),
153 fNPtBins(source.fNPtBins),
154 fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
155 fIsSelectedCandidate(source.fIsSelectedCandidate),
156 fSys(source.fSys),
157 fIsRejectSDDClusters(source.fIsRejectSDDClusters),
158 fFillGlobal(source.fFillGlobal)
159{
160 // Copy constructor
161}
162
163//________________________________________________________________________
164AliAnalysisTaskSED0Correlations::~AliAnalysisTaskSED0Correlations()
165{
166 if (fOutputMass) {
167 delete fOutputMass;
168 fOutputMass = 0;
169 }
170 if (fOutputCorr) {
171 delete fOutputCorr;
172 fOutputCorr = 0;
173 }
174 if (fOutputStudy) {
175 delete fOutputStudy;
176 fOutputStudy = 0;
177 }
178 if (fCutsD0) {
179 delete fCutsD0;
180 fCutsD0 = 0;
181 }
182 if (fNentries){
183 delete fNentries;
184 fNentries = 0;
185 }
186 if(fCounter){
187 delete fCounter;
188 fCounter=0;
189 }
190}
191
192//______________________________________________________________________________
193AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(const AliAnalysisTaskSED0Correlations& orig)
194{
195// Assignment
196 if (&orig == this) return *this; //if address is the same (same object), returns itself
197
198 AliAnalysisTaskSE::operator=(orig); //Uses the AliAnalysisTaskSE operator to assign the inherited part of the class
199 fNPtBinsCorr = orig.fNPtBinsCorr;
200 fBinLimsCorr = orig.fBinLimsCorr;
201 fPtThreshLow = orig.fPtThreshLow;
202 fPtThreshUp = orig.fPtThreshUp;
203 fEvents = orig.fEvents;
204 fAlreadyFilled = orig.fAlreadyFilled;
205 fOutputMass = orig.fOutputMass;
206 fOutputCorr = orig.fOutputCorr;
207 fOutputStudy = orig.fOutputStudy;
208 fNentries = orig.fNentries;
209 fCutsD0 = orig.fCutsD0;
210 fCutsTracks = orig.fCutsTracks;
211 fReadMC = orig.fReadMC;
212 fCounter = orig.fCounter;
213 fNPtBins = orig.fNPtBins;
214 fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
215 fIsSelectedCandidate = orig.fIsSelectedCandidate;
216 fSys = orig.fSys;
217 fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
218 fFillGlobal = orig.fFillGlobal;
219
220 return *this; //returns pointer of the class
221}
222
223//________________________________________________________________________
224void AliAnalysisTaskSED0Correlations::Init()
225{
226 // Initialization
227
228 if(fDebug > 1) printf("AnalysisTaskSED0Correlations::Init() \n");
229
230 //Copy of cuts objects
231 AliRDHFCutsD0toKpi* copyfCutsD0 = new AliRDHFCutsD0toKpi(*fCutsD0);
232 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
233 copyfCutsD0->SetName(nameoutput);
234
235 // Post the data
236 PostData(3,copyfCutsD0);
237 PostData(7,fCutsTracks);
238
239 return;
240}
241
242//________________________________________________________________________
243void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
244{
245
246 // Create the output container
247 //
248 if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
249
250 // Several histograms are more conveniently managed in a TList
251 fOutputMass = new TList();
252 fOutputMass->SetOwner();
253 fOutputMass->SetName("listMass");
254
255 fOutputCorr = new TList();
256 fOutputCorr->SetOwner();
257 fOutputCorr->SetName("correlationslist");
258
259 fOutputStudy = new TList();
260 fOutputStudy->SetOwner();
261 fOutputStudy->SetName("MCstudyplots");
262
263 TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ";
264
265 for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
266
267 nameMass="histMass_";
268 nameMass+=i;
269 nameSgn="histSgn_";
270 nameSgn+=i;
271 nameBkg="histBkg_";
272 nameBkg+=i;
273 nameRfl="histRfl_";
274 nameRfl+=i;
275
276 //histograms of invariant mass distributions
277
278 //MC signal
279 if(fReadMC){
280 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
281 tmpSt->Sumw2();
282
283 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
284 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
285 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
286 tmpBt->Sumw2();
287 tmpRt->Sumw2();
288 fOutputMass->Add(tmpSt);
289 fOutputMass->Add(tmpRt);
290 fOutputMass->Add(tmpBt);
291 }
292
293 //mass
294 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",120,1.5648,2.1648);
295 tmpMt->Sumw2();
296 fOutputMass->Add(tmpMt);
297 }
298
299 const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
300
301 fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex *** Integral(5,6) = pt out of bounds", 18,-0.5,17.5);
302
303 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
304 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
305 fReadMC ? fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected") : fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
306 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
307 fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
308 fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
309 if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
310 if(fReadMC && fSys==0){
311 fNentries->GetXaxis()->SetBinLabel(12,"K");
312 fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
313 }
314 fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
315 fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
316 if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
317 if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
318 fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
319 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
320
321 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
322 fCounter->Init();
323
324 CreateCorrelationsObjs(); //creates histos for correlations analysis
325
326 // Post the data
327 PostData(1,fOutputMass);
328 PostData(2,fNentries);
329 PostData(4,fCounter);
330 PostData(5,fOutputCorr);
331 PostData(6,fOutputStudy);
332
333 return;
334}
335
336//________________________________________________________________________
337void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
338{
339 // Execute analysis for current event:
340 // heavy flavor candidates association to MC truth
341 //cout<<"I'm in UserExec"<<endl;
342
343
344 //cuts order
345 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
346 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
347 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
348 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
349 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
350 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
351 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
352 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
353 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
354
355
356 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
357 fEvents++;
358
359 TString bname="D0toKpi";
360
361 TClonesArray *inputArray=0;
362
363 if(!aod && AODEvent() && IsStandardAOD()) {
364 // In case there is an AOD handler writing a standard AOD, use the AOD
365 // event in memory rather than the input (ESD) event.
366 aod = dynamic_cast<AliAODEvent*> (AODEvent());
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
372 if(aodHandler->GetExtensions()) {
373 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
374 AliAODEvent* aodFromExt = ext->GetAOD();
375 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
376 }
377 } else if(aod) {
378 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
379 }
380
381 if(!inputArray || !aod) {
382 printf("AliAnalysisTaskSED0Correlations::UserExec: input branch not found!\n");
383 return;
384 }
385
386 // fix for temporary bug in ESDfilter
387 // the AODs with null vertex pointer didn't pass the PhysSel
388 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
389
390 TClonesArray *mcArray = 0;
391 AliAODMCHeader *mcHeader = 0;
392
393 if(fReadMC) {
394 // load MC particles
395 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
396 if(!mcArray) {
397 printf("AliAnalysisTaskSED0Correlations::UserExec: MC particles branch not found!\n");
398 return;
399 }
400
401 // load MC header
402 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
403 if(!mcHeader) {
404 printf("AliAnalysisTaskSED0Correlations::UserExec: MC header branch not found!\n");
405 return;
406 }
407 }
408
409 //histogram filled with 1 for every AOD
410 fNentries->Fill(0);
411 fCounter->StoreEvent(aod,fCutsD0,fReadMC);
412
413 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
414 TString trigclass=aod->GetFiredTriggerClasses();
415 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
416
417 if(!fCutsD0->IsEventSelected(aod)) {
418 if(fCutsD0->GetWhyRejection()==1) // rejected for pileup
419 fNentries->Fill(13);
420 if(fSys==1 && (fCutsD0->GetWhyRejection()==2 || fCutsD0->GetWhyRejection()==3)) fNentries->Fill(15);
421 if(fCutsD0->GetWhyRejection()==7) fNentries->Fill(17);
422 return;
423 }
424
425 // Check the Nb of SDD clusters
426 if (fIsRejectSDDClusters) {
427 Bool_t skipEvent = kFALSE;
428 Int_t ntracks = 0;
429 if (aod) ntracks = aod->GetNTracks();
430 for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
431 // ... get the track
432 AliAODTrack * track = aod->GetTrack(itrack);
433 if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
434 skipEvent=kTRUE;
435 fNentries->Fill(16);
436 break;
437 }
438 }
439 if (skipEvent) return;
440 }
441
442 // AOD primary vertex
443 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
444
445 Bool_t isGoodVtx=kFALSE;
446
447 //vtx1->Print();
448 TString primTitle = vtx1->GetTitle();
449 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
450 isGoodVtx=kTRUE;
451 fNentries->Fill(3);
452 }
453
454 //Reset flag for tracks distributions fill
455 fAlreadyFilled=kFALSE;
456
457 // loop over candidates
458 Int_t nInD0toKpi = inputArray->GetEntriesFast();
459 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
460
461 if(fFillGlobal) { //loop on V0 and tracks for each event, to fill Pt distr. and InvMass distr.
462
463 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
464 Int_t pdgCodes[2] = {211,211};
465 Int_t idArrayV0[v0array->GetEntriesFast()][2];
466 for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
467 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
468 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
469 if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
470 if(fReadMC && (v0->MatchToMC(311,mcArray,2,pdgCodes)<0)) continue;
471 ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
472 ((TH1F*)fOutputStudy->FindObject("hist_Pt_K_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
473 }
474 }
475
476 for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tacks
477 AliAODTrack * track = aod->GetTrack(itrack);
478 //rejection of tracks
479 if(track->GetID() < 0) continue; //discard negative ID tracks
480 if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
481 if(!fCutsTracks->IsHadronSelected(track,vtx1,aod->GetMagneticField())) continue;
482 //pT distribution (in all events), charg and hadr cases
483 ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt());
484 if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
485 }
486
487 } //end of loops for global plot fill
488
489 Int_t nSelectedloose=0,nSelectedtight=0;
490 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
491 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
492
493 if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
494 fNentries->Fill(2);
495 continue; //skip the D0 from Dstar
496 }
497
498 if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
499 nSelectedloose++;
500 nSelectedtight++;
501 if(fSys==0){
502 if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
503 }
504 Int_t ptbin=fCutsD0->PtBin(d->Pt());
505 if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
506
507 fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
508 if(!fIsSelectedCandidate) continue;
509
510 if(!fReadMC) CalculateCorrelations(aod,d); //correlations on real data
511 else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
512 Int_t pdgDgD0toKpi[2]={321,211};
513 Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
514 if (labD0>-1) CalculateCorrelations(aod,d,labD0,mcArray);
515 }
516
517 FillMassHists(d,mcArray,fCutsD0,fOutputMass);
518 }
519
520 } //end for prongs
521 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
522 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
523
524 // Post the data
525 PostData(1,fOutputMass);
526 PostData(2,fNentries);
527 PostData(4,fCounter);
528 PostData(5,fOutputCorr);
529 PostData(6,fOutputStudy);
530
531 return;
532}
533
534//____________________________________________________________________________
535void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
536 //
537 // function used in UserExec to fill mass histograms:
538 //
539 if (!fIsSelectedCandidate) return;
540
541 if(fDebug>2) cout<<"Candidate selected"<<endl;
542
543 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
544 Int_t ptbin = cuts->PtBin(part->Pt());
545
546 TString fillthis="";
547 Int_t pdgDgD0toKpi[2]={321,211};
548 Int_t labD0=-1;
549 if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
550
551 //count candidates selected by cuts
552 fNentries->Fill(1);
553 //count true D0 selected by cuts
554 if (fReadMC && labD0>=0) fNentries->Fill(2);
555
556 if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
557
558 if(fReadMC){
559 if(labD0>=0) {
560 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
561 Int_t pdgD0 = partD0->GetPdgCode();
562 if (pdgD0==421){ //D0
563 fillthis="histSgn_";
564 fillthis+=ptbin;
565 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
566 } else{ //it was a D0bar
567 fillthis="histRfl_";
568 fillthis+=ptbin;
569 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
570 }
571 } else {//background
572 fillthis="histBkg_";
573 fillthis+=ptbin;
574 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
575 }
576 }else{
577 fillthis="histMass_";
578 fillthis+=ptbin;
579 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
580 }
581
582 }
583 if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
584
585 if(fReadMC){
586 if(labD0>=0) {
587 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
588 Int_t pdgD0 = partD0->GetPdgCode();
589
590 if (pdgD0==-421){ //D0bar
591 fillthis="histSgn_";
592 fillthis+=ptbin;
593 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
594 } else{
595 fillthis="histRfl_";
596 fillthis+=ptbin;
597 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
598 }
599 } else {//background or LS
600 fillthis="histBkg_";
601 fillthis+=ptbin;
602 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
603 }
604 }else{
605 fillthis="histMass_";
606 fillthis+=ptbin;
607 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
608 }
609 }
610
611 return;
612}
613
614//________________________________________________________________________
615void AliAnalysisTaskSED0Correlations::Terminate(Option_t */*option*/)
616{
617 // Terminate analysis
618 //
619 if(fDebug > 1) printf("AnalysisTaskSED0Correlations: Terminate() \n");
620
621 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
622 if (!fOutputMass) {
623 printf("ERROR: fOutputMass not available\n");
624 return;
625 }
626
627 fNentries = dynamic_cast<TH1F*>(GetOutputData(2));
628
629 if(!fNentries){
630 printf("ERROR: fNEntries not available\n");
631 return;
632 }
633
634 fCutsD0 = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(3));
635 if(!fCutsD0){
636 printf("ERROR: fCuts not available\n");
637 return;
638 }
639
640 fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(4));
641 if (!fCounter) {
642 printf("ERROR: fCounter not available\n");
643 return;
644 }
645 fOutputCorr = dynamic_cast<TList*> (GetOutputData(5));
646 if (!fOutputCorr) {
647 printf("ERROR: fOutputCorr not available\n");
648 return;
649 }
650 fOutputStudy = dynamic_cast<TList*> (GetOutputData(6));
651 if (!fOutputStudy) {
652 printf("ERROR: fOutputStudy not available\n");
653 return;
654 }
655 fCutsTracks = dynamic_cast<AliHFAssociatedTrackCuts*>(GetOutputData(7));
656 if(!fCutsTracks){
657 printf("ERROR: fCutsTracks not available\n");
658 return;
659 }
660
661 return;
662}
663
664//_________________________________________________________________________________________________
665Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
666 //
667 // checking whether the mother of the particles come from a charm or a bottom quark
668 //
669 printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
670
671 Int_t pdgGranma = 0;
672 Int_t mother = 0;
673 mother = mcPartCandidate->GetMother();
674 Int_t abspdgGranma =0;
675 Bool_t isFromB=kFALSE;
676 Bool_t isQuarkFound=kFALSE;
677
678 while (mother > 0){
679 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
680 if (mcGranma){
681 pdgGranma = mcGranma->GetPdgCode();
682 abspdgGranma = TMath::Abs(pdgGranma);
683 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
684 isFromB=kTRUE;
685 }
686 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
687 mother = mcGranma->GetMother();
688 }else{
689 AliError("Failed casting the mother particle!");
690 break;
691 }
692 }
693
694 if(isQuarkFound) {
695 if(isFromB) return 5;
696 else return 4;
697 }
698 else return 1;
699}
700
701//________________________________________________________________________
702void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
703//
704
705 TString namePlot = "";
706
707 //These for limits in THnSparse (one per bin, same limits).
708 //Vars: DeltaPhi, InvMass, PtTrack
709 Int_t nBinsPhi[4] = {32,150,30,3};
710 Double_t binMinPhi[4] = {-1.6,1.6,0.,0.}; //is the minimum for all the bins
711 Double_t binMaxPhi[4] = {4.8,2.2,3.0,3.}; //is the maximum for all the bins
712
713 for(Int_t i=0;i<fNPtBinsCorr;i++){
714
715 //THnSparse plots: correlations for various invariant mass (MC and data)
716 namePlot="hPhi_K_Bin";
717 namePlot+=i;
718
719 THnSparseI *hPhiK = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
720 hPhiK->Sumw2();
721 fOutputCorr->Add(hPhiK);
722
723 namePlot="hPhi_Kcharg_Bin";
724 namePlot+=i;
725
726 THnSparseI *hPhiH = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
727 hPhiH->Sumw2();
728 fOutputCorr->Add(hPhiH);
729
730 namePlot="hPhi_Charg_Bin";
731 namePlot+=i;
732
733 THnSparseI *hPhiC = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
734 hPhiC->Sumw2();
735 fOutputCorr->Add(hPhiC);
736
737 //histos for c/b origin for D0 (MC only)
738 if (fReadMC) {
739
740 //generic origin for tracks
741 namePlot="hPhi_K_From_c_Bin";
742 namePlot+=i;
743
744 THnSparseI *hPhiK_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
745 hPhiK_c->Sumw2();
746 fOutputCorr->Add(hPhiK_c);
747
748 namePlot="hPhi_Kcharg_From_c_Bin";
749 namePlot+=i;
750
751 THnSparseI *hPhiH_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
752 hPhiH_c->Sumw2();
753 fOutputCorr->Add(hPhiH_c);
754
755 namePlot="hPhi_Charg_From_c_Bin";
756 namePlot+=i;
757
758 THnSparseI *hPhiC_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
759 hPhiC_c->Sumw2();
760 fOutputCorr->Add(hPhiC_c);
761
762 namePlot="hPhi_K_From_b_Bin";
763 namePlot+=i;
764
765 THnSparseI *hPhiK_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
766 hPhiK_b->Sumw2();
767 fOutputCorr->Add(hPhiK_b);
768
769 namePlot="hPhi_Kcharg_From_b_Bin";
770 namePlot+=i;
771
772 THnSparseI *hPhiH_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
773 hPhiH_b->Sumw2();
774 fOutputCorr->Add(hPhiH_b);
775
776 namePlot="hPhi_Charg_From_b_Bin";
777 namePlot+=i;
778
779 THnSparseI *hPhiC_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
780 hPhiC_b->Sumw2();
781 fOutputCorr->Add(hPhiC_b);
782
783 //HF-only tracks (c for c->D0, b for b->D0)
784 namePlot="hPhi_K_HF_From_c_Bin";
785 namePlot+=i;
786
787 THnSparseI *hPhiK_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
788 hPhiK_HF_c->Sumw2();
789 fOutputCorr->Add(hPhiK_HF_c);
790
791 namePlot="hPhi_Kcharg_HF_From_c_Bin";
792 namePlot+=i;
793
794 THnSparseI *hPhiH_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
795 hPhiH_HF_c->Sumw2();
796 fOutputCorr->Add(hPhiH_HF_c);
797
798 namePlot="hPhi_Charg_HF_From_c_Bin";
799 namePlot+=i;
800 THnSparseI *hPhiC_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
801 hPhiC_HF_c->Sumw2();
802 fOutputCorr->Add(hPhiC_HF_c);
803
804 namePlot="hPhi_K_HF_From_b_Bin";
805 namePlot+=i;
806
807 THnSparseI *hPhiK_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
808 hPhiK_HF_b->Sumw2();
809 fOutputCorr->Add(hPhiK_HF_b);
810
811 namePlot="hPhi_Kcharg_HF_From_b_Bin";
812 namePlot+=i;
813
814 THnSparseI *hPhiH_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
815 hPhiH_HF_b->Sumw2();
816 fOutputCorr->Add(hPhiH_HF_b);
817
818 namePlot="hPhi_Charg_HF_From_b_Bin";
819 namePlot+=i;
820
821 THnSparseI *hPhiC_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
822 hPhiC_HF_b->Sumw2();
823 fOutputCorr->Add(hPhiC_HF_b);
824 }
825
826 //leading hadron correlations
827 namePlot="hPhi_Lead_Bin";
828 namePlot+=i;
829
830 TH2F *hCorrLead = new TH2F(namePlot.Data(), "Leading particle correlation; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
831 hCorrLead->Sumw2();
832 fOutputCorr->Add(hCorrLead);
833
834 if (fReadMC) {
835 namePlot="hPhi_Lead_From_c_Bin";
836 namePlot+=i;
837
838 TH2F *hCorrLead_c = new TH2F(namePlot.Data(), "Leading particle correlation - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
839 hCorrLead_c->Sumw2();
840 fOutputCorr->Add(hCorrLead_c);
841
842 namePlot="hPhi_Lead_From_b_Bin";
843 namePlot+=i;
844
845 TH2F *hCorrLead_b = new TH2F(namePlot.Data(), "Leading particle correlation - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
846 hCorrLead_b->Sumw2();
847 fOutputCorr->Add(hCorrLead_b);
848
849 namePlot="hPhi_Lead_HF_From_c_Bin";
850 namePlot+=i;
851
852 TH2F *hCorrLead_HF_c = new TH2F(namePlot.Data(), "Leading particle correlation HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
853 hCorrLead_HF_c->Sumw2();
854 fOutputCorr->Add(hCorrLead_HF_c);
855
856 namePlot="hPhi_Lead_HF_From_b_Bin";
857 namePlot+=i;
858
859 TH2F *hCorrLead_HF_b = new TH2F(namePlot.Data(), "Leading particle correlation HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
860 hCorrLead_HF_b->Sumw2();
861 fOutputCorr->Add(hCorrLead_HF_b);
862 }
863
864 //pT weighted correlations
865 namePlot="hPhi_Weig_Bin";
866 namePlot+=i;
867
868 TH2F *hCorrWeig = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
869 fOutputCorr->Add(hCorrWeig);
870
871 if (fReadMC) {
872 namePlot="hPhi_Weig_From_c_Bin";
873 namePlot+=i;
874
875 TH2F *hCorrWeig_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
876 fOutputCorr->Add(hCorrWeig_c);
877
878 namePlot="hPhi_Weig_From_b_Bin";
879 namePlot+=i;
880
881 TH2F *hCorrWeig_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
882 fOutputCorr->Add(hCorrWeig_b);
883
884 namePlot="hPhi_Weig_HF_From_c_Bin";
885 namePlot+=i;
886
887 TH2F *hCorrWeig_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
888 fOutputCorr->Add(hCorrWeig_HF_c);
889
890 namePlot="hPhi_Weig_HF_From_b_Bin";
891 namePlot+=i;
892
893 TH2F *hCorrWeig_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
894 fOutputCorr->Add(hCorrWeig_HF_b);
895 }
896
897 //pT weighted correlations (squared weights)
898 namePlot="hPhi_WeigSqr_Bin";
899 namePlot+=i;
900
901 TH2F *hCorrWeigSqr = new TH2F(namePlot.Data(), "Charged particle correlation (pT Weighted - squared); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
902 fOutputCorr->Add(hCorrWeigSqr);
903
904 if (fReadMC) {
905 namePlot="hPhi_WeigSqr_From_c_Bin";
906 namePlot+=i;
907
908 TH2F *hCorrWeigSqr_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
909 fOutputCorr->Add(hCorrWeigSqr_c);
910
911 namePlot="hPhi_WeigSqr_From_b_Bin";
912 namePlot+=i;
913
914 TH2F *hCorrWeigSqr_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
915 fOutputCorr->Add(hCorrWeigSqr_b);
916
917 namePlot="hPhi_WeigSqr_HF_From_c_Bin";
918 namePlot+=i;
919
920 TH2F *hCorrWeigSqr_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
921 fOutputCorr->Add(hCorrWeigSqr_HF_c);
922
923 namePlot="hPhi_WeigSqr_HF_From_b_Bin";
924 namePlot+=i;
925
926 TH2F *hCorrWeigSqr_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
927 fOutputCorr->Add(hCorrWeigSqr_HF_b);
928 }
929
930 //pT weighted correlations (inverse of pT distribution weights)
931 namePlot="hPhi_WeigDist_Bin";
932 namePlot+=i;
933
934 TH2F *hCorrWeigDist = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
935 fOutputCorr->Add(hCorrWeigDist);
936
937 if (fReadMC) {
938 namePlot="hPhi_WeigDist_From_c_Bin";
939 namePlot+=i;
940
941 TH2F *hCorrWeigDist_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
942 fOutputCorr->Add(hCorrWeigDist_c);
943
944 namePlot="hPhi_WeigDist_From_b_Bin";
945 namePlot+=i;
946
947 TH2F *hCorrWeigDist_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
948 fOutputCorr->Add(hCorrWeigDist_b);
949
950 namePlot="hPhi_WeigDist_HF_From_c_Bin";
951 namePlot+=i;
952
953 TH2F *hCorrWeigDist_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
954 fOutputCorr->Add(hCorrWeigDist_HF_c);
955
956 namePlot="hPhi_WeigDist_HF_From_b_Bin";
957 namePlot+=i;
958
959 TH2F *hCorrWeigDist_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
960 fOutputCorr->Add(hCorrWeigDist_HF_b);
961 }
962
963 //counters
964 namePlot = "hist_Count_Charg_Bin"; namePlot+=i;
965 TH1F *hCountC = new TH1F(namePlot.Data(), "Charged track counter; # Tracks",100,0.,100.);
966 hCountC->SetMinimum(0);
967 fOutputStudy->Add(hCountC);
968
969 namePlot = "hist_Count_Kcharg_Bin"; namePlot+=i;
970 TH1F *hCountH = new TH1F(namePlot.Data(), "Hadrons counter; # Tracks",100,0.,100.);
971 hCountH->SetMinimum(0);
972 fOutputStudy->Add(hCountH);
973
974 namePlot = "hist_Count_K_Bin"; namePlot+=i;
975 TH1F *hCountK = new TH1F(namePlot.Data(), "Kaons counter; # Tracks",10,0.,10.);
976 hCountK->SetMinimum(0);
977 fOutputStudy->Add(hCountK);
978
979 //pT distribution histos
980 namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
981 TH1F *hPtC = new TH1F(namePlot.Data(), "Charged track pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
982 hPtC->SetMinimum(0);
983 fOutputStudy->Add(hPtC);
984
985 namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
986 TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
987 hPtH->SetMinimum(0);
988 fOutputStudy->Add(hPtH);
989
990 namePlot = "hist_Pt_K_Bin"; namePlot+=i;
991 TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
992 hPtK->SetMinimum(0);
993 fOutputStudy->Add(hPtK);
994
995 //D* feeddown pions rejection histos
996 namePlot = "hDstarPions_Bin"; namePlot+=i;
997 TH1F *hDstarPions = new TH1F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.);
998 hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
999 hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
1000 hDstarPions->SetMinimum(0);
1001 fOutputStudy->Add(hDstarPions);
1002 }
1003
1004 //out of bin loop
1005 TH1F *hRejTracks = new TH1F("hRejTracks", "Tracks accepted/rejected; # Tracks",2,0.,2.);
1006 hRejTracks->SetMinimum(0);
1007 hRejTracks->GetXaxis()->SetBinLabel(1,"Accepted");
1008 hRejTracks->GetXaxis()->SetBinLabel(2,"Rejected");
1009 fOutputStudy->Add(hRejTracks);
1010
1011 if (fFillGlobal) { //all-events plots
1012 //pt distributions
1013 TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1014 hPtCAll->SetMinimum(0);
1015 fOutputStudy->Add(hPtCAll);
1016
1017 TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1018 hPtHAll->SetMinimum(0);
1019 fOutputStudy->Add(hPtHAll);
1020
1021 TH1F *hPtKAll = new TH1F("hist_Pt_K_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1022 hPtKAll->SetMinimum(0);
1023 fOutputStudy->Add(hPtKAll);
1024
1025 //K0 Invariant Mass plots
1026 TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1027 hK0MassInv->SetMinimum(0);
1028 fOutputStudy->Add(hK0MassInv);
1029 }
1030
1031 //for MC analysis only
1032 if (fReadMC) {
1033
1034 for(Int_t i=0;i<fNPtBinsCorr;i++){
1035
1036 //displacement histos
1037 namePlot="histDispl_K_Bin"; namePlot+=i;
1038 TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1039 hDisplK->SetMinimum(0);
1040 fOutputStudy->Add(hDisplK);
1041
1042 namePlot="histDispl_K_HF_Bin"; namePlot+=i;
1043 TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1044 hDisplK_HF->SetMinimum(0);
1045 fOutputStudy->Add(hDisplK_HF);
1046
1047 namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1048 TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1049 hDisplHadr->SetMinimum(0);
1050 fOutputStudy->Add(hDisplHadr);
1051
1052 namePlot="histDispl_Kcharg_HF_Bin"; namePlot+=i;
1053 TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1054 hDisplHadr_HF->SetMinimum(0);
1055 fOutputStudy->Add(hDisplHadr_HF);
1056
1057 namePlot="histDispl_Charg_Bin"; namePlot+=i;
1058 TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1059 hDisplCharg->SetMinimum(0);
1060 fOutputStudy->Add(hDisplCharg);
1061
1062 namePlot="histDispl_Charg_HF_Bin"; namePlot+=i;
1063 TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1064 hDisplCharg_HF->SetMinimum(0);
1065 fOutputStudy->Add(hDisplCharg_HF);
1066
1067 namePlot="histDispl_K_From_c_Bin"; namePlot+=i;
1068 TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1069 hDisplK_c->SetMinimum(0);
1070 fOutputStudy->Add(hDisplK_c);
1071
1072 namePlot="histDispl_K_HF_From_c_Bin"; namePlot+=i;
1073 TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1074 hDisplK_HF_c->SetMinimum(0);
1075 fOutputStudy->Add(hDisplK_HF_c);
1076
1077 namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1078 TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1079 hDisplHadr_c->SetMinimum(0);
1080 fOutputStudy->Add(hDisplHadr_c);
1081
1082 namePlot="histDispl_Kcharg_HF_From_c_Bin"; namePlot+=i;
1083 TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1084 hDisplHadr_HF_c->SetMinimum(0);
1085 fOutputStudy->Add(hDisplHadr_HF_c);
1086
1087 namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1088 TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1089 hDisplCharg_c->Sumw2();
1090 hDisplCharg_c->SetMinimum(0);
1091 fOutputStudy->Add(hDisplCharg_c);
1092
1093 namePlot="histDispl_Charg_HF_From_c_Bin"; namePlot+=i;
1094 TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1095 hDisplCharg_HF_c->SetMinimum(0);
1096 fOutputStudy->Add(hDisplCharg_HF_c);
1097
1098 namePlot="histDispl_K_From_b_Bin"; namePlot+=i;
1099 TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1100 hDisplK_b->SetMinimum(0);
1101 fOutputStudy->Add(hDisplK_b);
1102
1103 namePlot="histDispl_K_HF_From_b_Bin"; namePlot+=i;
1104 TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1105 hDisplK_HF_b->SetMinimum(0);
1106 fOutputStudy->Add(hDisplK_HF_b);
1107
1108 namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1109 TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1110 hDisplHadr_b->SetMinimum(0);
1111 fOutputStudy->Add(hDisplHadr_b);
1112
1113 namePlot="histDispl_Kcharg_HF_From_b_Bin"; namePlot+=i;
1114 TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1115 hDisplHadr_HF_b->SetMinimum(0);
1116 fOutputStudy->Add(hDisplHadr_HF_b);
1117
1118 namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1119 TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1120 hDisplCharg_b->SetMinimum(0);
1121 fOutputStudy->Add(hDisplCharg_b);
1122
1123 namePlot="histDispl_Charg_HF_From_b_Bin"; namePlot+=i;
1124 TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1125 hDisplCharg_HF_b->SetMinimum(0);
1126 fOutputStudy->Add(hDisplCharg_HF_b);
1127
1128 //origin of tracks histos
1129 namePlot="histOrig_Charg_Bin"; namePlot+=i;
1130 TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1131 hOrigin_Charm->SetMinimum(0);
1132 hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1133 hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1134 hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1135 hOrigin_Charm->GetXaxis()->SetBinLabel(4,"B->#");
1136 hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1137 hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->D->#");
1138 hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1139 hOrigin_Charm->GetXaxis()->SetBinLabel(8,"c hadr.");
1140 hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1141 fOutputStudy->Add(hOrigin_Charm);
1142
1143 namePlot="histOrig_Kcharg_Bin"; namePlot+=i;
1144 TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1145 hOrigin_Kcharg->SetMinimum(0);
1146 hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1147 hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1148 hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1149 hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"B->#");
1150 hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1151 hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->D->#");
1152 hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1153 hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"c hadr.");
1154 hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1155 fOutputStudy->Add(hOrigin_Kcharg);
1156
1157 namePlot="histOrig_K_Bin"; namePlot+=i;
1158 TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1159 hOrigin_K->SetMinimum(0);
1160 hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1161 hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1162 hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1163 hOrigin_K->GetXaxis()->SetBinLabel(4,"B->#");
1164 hOrigin_K->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1165 hOrigin_K->GetXaxis()->SetBinLabel(6,"B->D->#");
1166 hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1167 hOrigin_K->GetXaxis()->SetBinLabel(8,"c hadr.");
1168 hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1169 fOutputStudy->Add(hOrigin_K);
1170
1171 //origin of D0 histos
1172 namePlot="histOrig_D0_Bin"; namePlot+=i;
1173 TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1174 hOrigin_D0->SetMinimum(0);
1175 hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1176 hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1177 fOutputStudy->Add(hOrigin_D0);
1178 }
1179 }
1180
1181}
1182
1183//________________________________________________________________________
1184void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODEvent* aod, AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1185//
1186// Method for correlations D0-hadrons study
1187//
1188
1189 Double_t mD0, mD0bar, deltaphi, d0, d0err;
1190 d->InvMassD0(mD0,mD0bar);
1191 Int_t ptbin = PtBinCorr(d->Pt());
1192 if(ptbin < 0) return;
1193 AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
1194
1195 Int_t N_Charg = 0, N_Kcharg = 0, N_Kaons = 0;
1196
1197 if(fReadMC == 0) {
1198 Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1199 Double_t highPt = 0; Double_t lead[2] = {0,0}; //infos for leading particle (pt,deltaphi)
1200
1201 for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
1202 AliAODTrack * track = aod->GetTrack(itrack);
1203
1204 if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
1205
1206 GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
1207
1208 ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
1209 deltaphi = d->Phi()-track->Phi();
1210 if (deltaphi < -1.571) deltaphi+=6.283;
1211 if (deltaphi > 4.712) deltaphi-=6.283;
1212 Double_t pttrack = track->Pt();
1213
1214 if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi;} //for leading particle
1215
1216 //Lines needed to include overflow into THnSparse projections!
1217 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1218 Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
1219 if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
1220 if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
1221
1222 //variables for filling histos
1223 Double_t fillSpPhiD0[4] = {deltaphi,mD0,pttrack,d0/d0err};
1224 Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,pttrack,d0/d0err};
1225
1226 //generic charged tracks (NO PID selection)
1227 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1228 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1229 }
1230 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1231 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1232 }
1233 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
1234 N_Charg++;
1235 //pT_weighted track correlations fill
1236 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1237 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
1238 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1239 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1240 }
1241 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1242 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1243 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1244 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1245 }
1246
1247 //hadron identification
1248 if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) {
1249 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1250 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1251 }
1252 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1253 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1254 }
1255 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
1256 N_Kcharg++;
1257 } //end hadron case
1258
1259 } //end tracks loop
1260
1261 //kaon identification
1262 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
1263 Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
1264 for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
1265 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
1266
1267 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
1268 if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1269 if(SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
1270
1271 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1272 Double_t ptV0=v0->Pt(), deltaphiV0=d->Phi()-v0->Phi();
1273 if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
1274 deltaphiV0 = d->Phi()-v0->Phi();
1275 if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
1276 if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
1277
1278 Double_t fillSpPhiD0K0[4] = {deltaphiV0,mD0,ptV0,0.};
1279 Double_t fillSpPhiD0barK0[4] = {deltaphiV0,mD0bar,ptV0,0.};
1280
1281 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1282 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
1283 }
1284 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1285 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
1286 }
1287 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K_Bin%d",ptbin)))->Fill(v0->Pt());
1288 N_Kaons++;
1289 } // end kaon case
1290
1291 //Leading track correlations fill
1292 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1293 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0);
1294 }
1295 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1296 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
1297 }
1298
1299 } //end fReadMC = 0
1300
1301 if(fReadMC == 1) {
1302 Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1303 Double_t highPt = 0; Double_t lead[3] = {0,0,0}; //infos for leading particle (pt,deltaphi,origtrack)
1304
1305 //Origin of D0
1306 TString orig=""; Int_t origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1307 switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0)))
1308 {
1309 case 4:
1310 orig = "_From_c";
1311 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1312 break;
1313 case 5:
1314 orig = "_From_b";
1315 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1316 break;
1317 default:
1318 return;
1319 }
1320
1321 for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
1322 AliAODTrack * track = aod->GetTrack(itrack);
1323
1324 if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
1325
1326 Int_t label = track->GetLabel();
1327 if (label<0) continue; //discard negative label tracks
1328 Int_t PDGtrack = ((AliAODMCParticle*)mcArray->At(label))->GetPdgCode();
1329
1330 GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
1331
1332 //Infos on track (origin, phi, eta)
1333 Int_t origTr = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(label));
1334 deltaphi = d->Phi()-track->Phi();
1335 if (deltaphi < -1.571) deltaphi+=6.283;
1336 if (deltaphi > 4.712) deltaphi-=6.283;
1337 Double_t pttrack = track->Pt();
1338
1339 if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi; lead[2] = origTr;} //for leading particle
1340
1341 //Lines needed to include overflow into THnSparse projections!
1342 Double_t d0orig = d0;
1343 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1344 Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
1345 if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
1346 if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
1347
1348 //variables for filling histos
1349 Double_t fillSpPhiD0[4] = {deltaphi,mD0,pttrack,d0/d0err};
1350 Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,pttrack,d0/d0err};
1351
1352 //generic charged tracks case
1353 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
1354 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1355 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1356 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//c tr
1357 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//b tr
1358 }
1359 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1360 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1361 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1362 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1363 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1364 }
1365 if(!fAlreadyFilled) {
1366 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
1367 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF_Bin%d",ptbin)))->Fill(d0orig);
1368 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
1369 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1370 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
1371 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Charg_Bin%d",ptbin)))->Fill(origTr);
1372 }
1373 N_Charg++;
1374 //pT_weighted track correlations fill
1375 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1376 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
1377 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);
1378 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1379 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1380 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1381 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1382 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);//c tr
1383 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);//b tr
1384 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));//c tr
1385 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));//b tr
1386 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));//c tr
1387 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));//b tr
1388 }
1389 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1390 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1391 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1392 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1393 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1394 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1395 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1396 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);//c tr
1397 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);//b tr
1398 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));//c tr
1399 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));//b tr
1400 if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));//c tr
1401 if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));//b tr
1402 }
1403
1404 //hadron identification (K/pi/p from MCTruth)
1405 if(TMath::Abs(PDGtrack) == 321) {
1406 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
1407 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1408 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1409 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1410 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1411 }
1412 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1413 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1414 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1415 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1416 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1417 }
1418 if(!fAlreadyFilled) {
1419 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
1420 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF_Bin%d",ptbin)))->Fill(d0orig);
1421 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
1422 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1423 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
1424 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Kcharg_Bin%d",ptbin)))->Fill(origTr);
1425 }
1426 N_Kcharg++;
1427
1428 } //end hadron case
1429
1430 } //end tracks loop
1431
1432 //Kaon identification
1433 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
1434 Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
1435 for(int iV0=0; iV0<v0array->GetEntriesFast();iV0++) {
1436 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
1437 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
1438
1439 if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1440 if(SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
1441 Int_t pdgCodes[2] = {211,211};
1442 Int_t labV0 = v0->MatchToMC(311,mcArray,2,pdgCodes);
1443 if(labV0<=0) continue;
1444
1445 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1446 Double_t ptV0=v0->Pt(), deltaphiV0=d->Phi()-v0->Phi();
1447 deltaphiV0 = d->Phi()-v0->Phi();
1448 if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
1449 if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
1450 if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
1451
1452 Double_t fillSpPhiD0K0[4] = {deltaphiV0,mD0,ptV0,0.};
1453 Double_t fillSpPhiD0barK0[4] = {deltaphiV0,mD0bar,ptV0,0.};
1454
1455 Int_t origV0 = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(labV0));
1456
1457 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
1458 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
1459 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
1460 if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
1461 if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
1462 }
1463 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1464 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
1465 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
1466 if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
1467 if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
1468 }
1469 if(!fAlreadyFilled) {
1470 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_Bin%d",ptbin)))->Fill(0.); //Fills displacement histos
1471 if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_HF_Bin%d",ptbin)))->Fill(0.);
1472 if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(0.);
1473 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K%s_Bin%d",orig.Data(),ptbin)))->Fill(0.); //Fills displacement histos
1474 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K_Bin%d",ptbin)))->Fill(v0->Pt());
1475 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_K_Bin%d",ptbin)))->Fill(origV0);
1476 }
1477 N_Kaons++;
1478 } // end kaon case
1479
1480 //leading track correlations fill
1481 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1482 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0); //c and b D0
1483 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0); //c or b D0
1484 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0);
1485 if(origD0==5&&(int)lead[2]>=3&&(int)lead[2]<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0);
1486 }
1487 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1488 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
1489 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar); //c or b D0
1490 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar);
1491 if(origD0==5&&(int)lead[2]>=3&&(int)lead[2]<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar);
1492 }
1493
1494 } //end fReadMC = 1
1495
1496 if (!fAlreadyFilled) {
1497 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Charg_Bin%d",ptbin)))->Fill(N_Charg);
1498 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Kcharg_Bin%d",ptbin)))->Fill(N_Kcharg);
1499 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_K_Bin%d",ptbin)))->Fill(N_Kaons);
1500 }
1501
1502 fAlreadyFilled=kTRUE; //distribution plots for tracks filled
1503
1504}
1505
1506//_________________________________________________________________________________________________
1507Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1508 //
1509 // checks on particle (#) origin:
1510 // 0) Not HF
1511 // 1) D->#
1512 // 2) D->X->#
1513 // 3) B->#
1514 // 4) B->X-># (X!=D)
1515 // 5) B->D->#
1516 // 6) B->D->X->#
1517 // 7) c hadronization
1518 // 8) b hadronization
1519 //
1520 printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
1521
1522 Int_t pdgGranma = 0;
1523 Int_t mother = 0;
1524 mother = mcPartCandidate->GetMother();
1525 Int_t istep = 0;
1526 Int_t abspdgGranma =0;
1527 Bool_t isFromB=kFALSE;
1528 Bool_t isDdaugh=kFALSE;
1529 Bool_t isDchaindaugh=kFALSE;
1530 Bool_t isBdaugh=kFALSE;
1531 Bool_t isBchaindaugh=kFALSE;
1532 Bool_t isQuarkFound=kFALSE;
1533
1534 while (mother > 0){
1535 istep++;
1536 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1537 if (mcGranma){
1538 pdgGranma = mcGranma->GetPdgCode();
1539 abspdgGranma = TMath::Abs(pdgGranma);
1540 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1541 isBchaindaugh=kTRUE;
1542 if(istep==1) isBdaugh=kTRUE;
1543 }
1544 if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
1545 isDchaindaugh=kTRUE;
1546 if(istep==1) isDdaugh=kTRUE;
1547 }
1548 if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
1549 mother = mcGranma->GetMother();
1550 }else{
1551 AliError("Failed casting the mother particle!");
1552 break;
1553 }
1554 }
1555
1556 //decides what to return based on the flag status
1557 if(isQuarkFound) {
1558 if(!isFromB) { //charm
1559 if(isDdaugh) return 1; //charm immediate
1560 else if(isDchaindaugh) return 2; //charm chain
1561 else return 7; //charm hadronization
1562 }
1563 else { //beauty
1564 if(isBdaugh) return 3; //b immediate
1565 else if(isBchaindaugh) { //b chain
1566 if(isDchaindaugh) {
1567 if(isDdaugh) return 5; //d immediate
1568 return 6; //d chain
1569 }
1570 else return 4; //b, not d
1571 }
1572 else return 8; //b hadronization
1573 }
1574 }
1575 else return 0; //no HF quark
1576}
1577
1578
1579//________________________________________________________________________
1580Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
1581 //
1582 //give the pt bin where the pt lies.
1583 //
1584 Int_t ptbin=-1;
1585 if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
1586
1587 Int_t i = 0;
1588 while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
1589
1590 return ptbin;
1591}
1592
1593//________________________________________________________________________
1594Double_t AliAnalysisTaskSED0Correlations::PtWeig(Int_t ptbin, Double_t x) const {
1595 //
1596 //gives the weight for Weighted Corrs (inverse of pT distribution, from 3 D0 pt bins)
1597 //
1598 if(x>11.5) x=11.5; if(x<0.3) x=0.3; //bounds for fit of distributions!
1599 if(ptbin >= 3 && ptbin <= 5) return 1/(0.22958*(1.507e+03-5.035e+02*x+5.681e+01*x*x-2.186e+00*x*x*x+exp(1.336e+01-2.146e+00*x+1.206e-01*x*x)));
1600 if(ptbin >= 6 && ptbin <= 8) return 1/(1.71828*(1.984e+02-6.113e+01*x+6.444e+00*x*x-2.342e-01*x*x*x+exp(1.023e+01-2.059e+00*x+1.194e-01*x*x)));
1601 if(ptbin >= 9 && ptbin <= 10) return 1/(1.25905*(1.990e+02-6.306e+01*x+6.995e+00*x*x-2.681e-01*x*x*x+exp(9.125e+00-2.053e+00*x+1.276e-01*x*x)));
1602
1603 return 0; //for other bins plot is disabled!
1604}
1605
1606//________________________________________________________________________
1607void AliAnalysisTaskSED0Correlations::GetImpParameter(AliAODTrack *track, AliAODEvent* aod, Double_t &d0, Double_t &d0err) const {
1608 //
1609 //calculates d0 and error on d0 for the track
1610 //
1611 Double_t d0z0[2],covd0z0[3];
1612 AliESDtrack esdTrack(track); // AliESDTrack conversion of AOD track
1613 esdTrack.PropagateToDCA(aod->GetPrimaryVertex(),aod->GetMagneticField(),10000.,d0z0,covd0z0);
1614 d0 = TMath::Abs(d0z0[0]); // impact parameter
1615 d0err = TMath::Sqrt(covd0z0[0]); // resolution of impact parameter
1616}
1617
1618//________________________________________________________________________
1619Bool_t AliAnalysisTaskSED0Correlations::TrackSelectedInLoop(AliAODTrack* track, AliAODRecoDecayHF2Prong *d, AliAODEvent *aod, Int_t ptbin, Int_t idDaughs[]) const {
1620 //
1621 //rejection of tracks in loop for correlations
1622 //
1623 Bool_t output = kTRUE;
1624 AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
1625 Double_t bz = aod->GetMagneticField();
1626 if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1] || track->GetID() < 0) output = kFALSE; //discards daughters of candidate
1627 if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) output = kFALSE; //discard tracks outside pt range for hadrons/K
1628 if(!fCutsTracks->IsHadronSelected(track,vtx,bz)) { //track discarded by quality cuts
1629 ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(1);
1630 output = kFALSE;
1631 } else ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
1632 if(!fCutsTracks->InvMassDstarRejection(d,track,fIsSelectedCandidate)) {
1633 ((TH1F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1);
1634 output = kFALSE;
1635 } else ((TH1F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0); //rejects possible pions from D* using inv.mass
1636
1637 return output;
1638}
1639
1640//---------------------------------------------------------------------------
1641Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
1642{
1643 //
1644 // Selection for K0 hypotheses
1645 // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
1646 // 2 = no previous selections
1647
1648 if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
1649
1650 AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
1651 AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
1652
1653 if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
1654 if(v0->Pt() < 3. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0040) return kFALSE;
1655 if(v0->Pt() > 3. && v0->Pt() < 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0052) return kFALSE;
1656 if(v0->Pt() > 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0075) return kFALSE;
1657 }
1658
1659 //This part removes double counting for swapped tracks!
1660 Int_t i = 0; //while loop (until the last-written entry pair of ID!
1661 while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
1662 if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
1663 (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
1664 i++;
1665 }
1666 idArrayV0[i][0]=v0Daug1->GetID();
1667 idArrayV0[i][1]=v0Daug2->GetID();
1668
1669 return kTRUE;
1670}
1671
1672//________________________________________________________________________
1673void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
1674
1675 cout << "--------------------------\n";
1676 cout << "PtBins = " << fNPtBinsCorr << "\n";
1677 cout << "PtBin limits--------------\n";
1678 for (int i=0; i<fNPtBinsCorr; i++) {
1679 cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
1680 }
1681 cout << "\n--------------------------\n";
1682 cout << "PtBin tresh. tracks low---\n";
1683 for (int i=0; i<fNPtBinsCorr; i++) {
1684 cout << fPtThreshLow.at(i) << "; ";
1685 }
1686 cout << "PtBin tresh. tracks up----\n";
1687 for (int i=0; i<fNPtBinsCorr; i++) {
1688 cout << fPtThreshUp.at(i) << "; ";
1689 }
1690 cout << "\n--------------------------\n";
1691 cout << "MC Truth = "<<fReadMC<<"\n";
1692 cout << "--------------------------\n";
1693}
1694