]> 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
increase class def
[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
f80e7bba 472 ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
a2ad7da1 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
f80e7bba 709 Int_t nBinsPhi[5] = {32,150,30,3,16};
710 Double_t binMinPhi[5] = {-1.6,1.6,0.,0.,-1.6}; //is the minimum for all the bins
711 Double_t binMaxPhi[5] = {4.8,2.2,3.0,3.,1.6}; //is the maximum for all the bins
a2ad7da1 712
713 for(Int_t i=0;i<fNPtBinsCorr;i++){
714
715 //THnSparse plots: correlations for various invariant mass (MC and data)
f80e7bba 716 namePlot="hPhi_K0_Bin";
a2ad7da1 717 namePlot+=i;
718
f80e7bba 719 THnSparseI *hPhiK = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 720 hPhiK->Sumw2();
721 fOutputCorr->Add(hPhiK);
722
723 namePlot="hPhi_Kcharg_Bin";
724 namePlot+=i;
725
f80e7bba 726 THnSparseI *hPhiH = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 727 hPhiH->Sumw2();
728 fOutputCorr->Add(hPhiH);
729
730 namePlot="hPhi_Charg_Bin";
731 namePlot+=i;
732
f80e7bba 733 THnSparseI *hPhiC = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 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
f80e7bba 741 namePlot="hPhi_K0_From_c_Bin";
a2ad7da1 742 namePlot+=i;
743
f80e7bba 744 THnSparseI *hPhiK_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 745 hPhiK_c->Sumw2();
746 fOutputCorr->Add(hPhiK_c);
747
748 namePlot="hPhi_Kcharg_From_c_Bin";
749 namePlot+=i;
750
f80e7bba 751 THnSparseI *hPhiH_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 752 hPhiH_c->Sumw2();
753 fOutputCorr->Add(hPhiH_c);
754
755 namePlot="hPhi_Charg_From_c_Bin";
756 namePlot+=i;
757
f80e7bba 758 THnSparseI *hPhiC_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 759 hPhiC_c->Sumw2();
760 fOutputCorr->Add(hPhiC_c);
761
f80e7bba 762 namePlot="hPhi_K0_From_b_Bin";
a2ad7da1 763 namePlot+=i;
764
f80e7bba 765 THnSparseI *hPhiK_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 766 hPhiK_b->Sumw2();
767 fOutputCorr->Add(hPhiK_b);
768
769 namePlot="hPhi_Kcharg_From_b_Bin";
770 namePlot+=i;
771
f80e7bba 772 THnSparseI *hPhiH_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 773 hPhiH_b->Sumw2();
774 fOutputCorr->Add(hPhiH_b);
775
776 namePlot="hPhi_Charg_From_b_Bin";
777 namePlot+=i;
778
f80e7bba 779 THnSparseI *hPhiC_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 780 hPhiC_b->Sumw2();
781 fOutputCorr->Add(hPhiC_b);
782
783 //HF-only tracks (c for c->D0, b for b->D0)
f80e7bba 784 namePlot="hPhi_K0_HF_From_c_Bin";
a2ad7da1 785 namePlot+=i;
786
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 788 hPhiK_HF_c->Sumw2();
789 fOutputCorr->Add(hPhiK_HF_c);
790
791 namePlot="hPhi_Kcharg_HF_From_c_Bin";
792 namePlot+=i;
793
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 795 hPhiH_HF_c->Sumw2();
796 fOutputCorr->Add(hPhiH_HF_c);
797
798 namePlot="hPhi_Charg_HF_From_c_Bin";
799 namePlot+=i;
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 801 hPhiC_HF_c->Sumw2();
802 fOutputCorr->Add(hPhiC_HF_c);
803
f80e7bba 804 namePlot="hPhi_K0_HF_From_b_Bin";
a2ad7da1 805 namePlot+=i;
806
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 808 hPhiK_HF_b->Sumw2();
809 fOutputCorr->Add(hPhiK_HF_b);
810
811 namePlot="hPhi_Kcharg_HF_From_b_Bin";
812 namePlot+=i;
813
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 815 hPhiH_HF_b->Sumw2();
816 fOutputCorr->Add(hPhiH_HF_b);
817
818 namePlot="hPhi_Charg_HF_From_b_Bin";
819 namePlot+=i;
820
f80e7bba 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)",5,nBinsPhi,binMinPhi,binMaxPhi);
a2ad7da1 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
f80e7bba 974 namePlot = "hist_Count_K0_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 990 namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
a2ad7da1 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;
f80e7bba 997 TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,300,1.6,2.2);
a2ad7da1 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
f80e7bba 1021 TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
a2ad7da1 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
f80e7bba 1037 namePlot="histDispl_K0_Bin"; namePlot+=i;
a2ad7da1 1038 TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1039 hDisplK->SetMinimum(0);
1040 fOutputStudy->Add(hDisplK);
1041
f80e7bba 1042 namePlot="histDispl_K0_HF_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1067 namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1072 namePlot="histDispl_K0_HF_From_c_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1098 namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1103 namePlot="histDispl_K0_HF_From_b_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1157 namePlot="histOrig_K0_Bin"; namePlot+=i;
a2ad7da1 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
f80e7bba 1189 Double_t mD0, mD0bar, deltaphi, d0, d0err, deltaeta;
a2ad7da1 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();
f80e7bba 1210 deltaeta = d->Eta()-track->Eta();
a2ad7da1 1211 if (deltaphi < -1.571) deltaphi+=6.283;
1212 if (deltaphi > 4.712) deltaphi-=6.283;
1213 Double_t pttrack = track->Pt();
1214
1215 if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi;} //for leading particle
1216
1217 //Lines needed to include overflow into THnSparse projections!
1218 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1219 Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
f80e7bba 1220 Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
a2ad7da1 1221 if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
1222 if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
f80e7bba 1223 if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1224 if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
a2ad7da1 1225
1226 //variables for filling histos
f80e7bba 1227 Double_t fillSpPhiD0[5] = {deltaphi,mD0,pttrack,d0/d0err,deltaeta};
1228 Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,pttrack,d0/d0err,deltaeta};
a2ad7da1 1229
1230 //generic charged tracks (NO PID selection)
1231 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1232 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1233 }
1234 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1235 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1236 }
1237 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
1238 N_Charg++;
1239 //pT_weighted track correlations fill
1240 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1241 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
1242 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1243 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1244 }
1245 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1246 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1247 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1248 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1249 }
1250
f80e7bba 1251 //Charged Kaon identification
a2ad7da1 1252 if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) {
1253 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1254 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1255 }
1256 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1257 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1258 }
1259 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
1260 N_Kcharg++;
1261 } //end hadron case
1262
1263 } //end tracks loop
1264
1265 //kaon identification
1266 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
1267 Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
1268 for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
1269 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
1270
1271 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
1272 if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
f80e7bba 1273 if(!SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
a2ad7da1 1274
f80e7bba 1275 Double_t ptV0=v0->Pt(), deltaphiV0, deltaetaV0;
1276 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1277 Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
a2ad7da1 1278 if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
1279 deltaphiV0 = d->Phi()-v0->Phi();
f80e7bba 1280 deltaetaV0 = d->Eta()-v0->Eta();
a2ad7da1 1281 if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
1282 if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
f80e7bba 1283 if(deltaetaV0 > EtaLim_Sparse) deltaetaV0 = EtaLim_Sparse-0.01;
1284 if(deltaetaV0 < -EtaLim_Sparse) deltaetaV0 = -EtaLim_Sparse+0.01;
a2ad7da1 1285
f80e7bba 1286 Double_t fillSpPhiD0K0[5] = {deltaphiV0,mD0,ptV0,0.,deltaetaV0};
1287 Double_t fillSpPhiD0barK0[5] = {deltaphiV0,mD0bar,ptV0,0.,deltaetaV0};
a2ad7da1 1288
1289 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
f80e7bba 1290 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
a2ad7da1 1291 }
1292 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
f80e7bba 1293 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
a2ad7da1 1294 }
f80e7bba 1295 if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K0_Bin%d",ptbin)))->Fill(v0->Pt());
a2ad7da1 1296 N_Kaons++;
1297 } // end kaon case
1298
1299 //Leading track correlations fill
1300 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1301 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0);
1302 }
1303 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1304 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
1305 }
1306
1307 } //end fReadMC = 0
1308
1309 if(fReadMC == 1) {
1310 Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1311 Double_t highPt = 0; Double_t lead[3] = {0,0,0}; //infos for leading particle (pt,deltaphi,origtrack)
1312
1313 //Origin of D0
1314 TString orig=""; Int_t origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1315 switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0)))
1316 {
1317 case 4:
1318 orig = "_From_c";
1319 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1320 break;
1321 case 5:
1322 orig = "_From_b";
1323 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1324 break;
1325 default:
1326 return;
1327 }
1328
1329 for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
1330 AliAODTrack * track = aod->GetTrack(itrack);
1331
1332 if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
1333
1334 Int_t label = track->GetLabel();
1335 if (label<0) continue; //discard negative label tracks
1336 Int_t PDGtrack = ((AliAODMCParticle*)mcArray->At(label))->GetPdgCode();
1337
1338 GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
1339
1340 //Infos on track (origin, phi, eta)
1341 Int_t origTr = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(label));
1342 deltaphi = d->Phi()-track->Phi();
f80e7bba 1343 deltaeta = d->Eta()-track->Eta();
a2ad7da1 1344 if (deltaphi < -1.571) deltaphi+=6.283;
1345 if (deltaphi > 4.712) deltaphi-=6.283;
1346 Double_t pttrack = track->Pt();
1347
1348 if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi; lead[2] = origTr;} //for leading particle
1349
1350 //Lines needed to include overflow into THnSparse projections!
1351 Double_t d0orig = d0;
1352 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1353 Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
f80e7bba 1354 Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
a2ad7da1 1355 if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
1356 if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
f80e7bba 1357 if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1358 if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
a2ad7da1 1359
1360 //variables for filling histos
f80e7bba 1361 Double_t fillSpPhiD0[5] = {deltaphi,mD0,pttrack,d0/d0err,deltaeta};
1362 Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,pttrack,d0/d0err,deltaeta};
a2ad7da1 1363
1364 //generic charged tracks case
1365 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
1366 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1367 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1368 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//c tr
1369 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//b tr
1370 }
1371 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1372 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1373 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1374 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1375 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1376 }
1377 if(!fAlreadyFilled) {
1378 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
1379 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF_Bin%d",ptbin)))->Fill(d0orig);
1380 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
1381 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1382 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
1383 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Charg_Bin%d",ptbin)))->Fill(origTr);
1384 }
1385 N_Charg++;
1386 //pT_weighted track correlations fill
1387 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1388 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
1389 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);
1390 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1391 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
1392 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1393 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
1394 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
1395 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
1396 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
1397 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
1398 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
1399 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
1400 }
1401 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1402 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1403 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);
1404 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1405 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
1406 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1407 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
1408 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
1409 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
1410 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
1411 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
1412 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
1413 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
1414 }
1415
f80e7bba 1416 //Charged Kaon identification (K from MCTruth)
a2ad7da1 1417 if(TMath::Abs(PDGtrack) == 321) {
1418 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
1419 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
1420 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1421 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1422 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1423 }
1424 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1425 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
1426 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1427 if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1428 if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1429 }
1430 if(!fAlreadyFilled) {
1431 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
1432 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF_Bin%d",ptbin)))->Fill(d0orig);
1433 if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
1434 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1435 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
1436 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Kcharg_Bin%d",ptbin)))->Fill(origTr);
1437 }
1438 N_Kcharg++;
1439
1440 } //end hadron case
1441
1442 } //end tracks loop
1443
1444 //Kaon identification
1445 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
1446 Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
1447 for(int iV0=0; iV0<v0array->GetEntriesFast();iV0++) {
1448 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
1449 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
1450
1451 if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
f80e7bba 1452 if(!SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
a2ad7da1 1453 Int_t pdgCodes[2] = {211,211};
1454 Int_t labV0 = v0->MatchToMC(311,mcArray,2,pdgCodes);
1455 if(labV0<=0) continue;
1456
f80e7bba 1457 Double_t ptV0=v0->Pt(), deltaphiV0, deltaetaV0;
1458 Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1459 Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
a2ad7da1 1460 deltaphiV0 = d->Phi()-v0->Phi();
f80e7bba 1461 deltaetaV0 = d->Eta()-v0->Eta();
a2ad7da1 1462 if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
1463 if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
1464 if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
f80e7bba 1465 if(deltaetaV0 > EtaLim_Sparse) deltaetaV0 = EtaLim_Sparse-0.01;
1466 if(deltaetaV0 < -EtaLim_Sparse) deltaetaV0 = -EtaLim_Sparse+0.01;
a2ad7da1 1467
f80e7bba 1468 Double_t fillSpPhiD0K0[5] = {deltaphiV0,mD0,ptV0,0.,deltaetaV0};
1469 Double_t fillSpPhiD0barK0[5] = {deltaphiV0,mD0bar,ptV0,0.,deltaetaV0};
a2ad7da1 1470
1471 Int_t origV0 = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(labV0));
1472
1473 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
f80e7bba 1474 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
1475 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
1476 if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
1477 if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
a2ad7da1 1478 }
1479 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
f80e7bba 1480 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
1481 ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
1482 if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
1483 if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
a2ad7da1 1484 }
1485 if(!fAlreadyFilled) {
f80e7bba 1486 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_Bin%d",ptbin)))->Fill(0.); //Fills displacement histos
1487 if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_HF_Bin%d",ptbin)))->Fill(0.);
1488 if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(0.);
1489 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(0.); //Fills displacement histos
1490 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K0_Bin%d",ptbin)))->Fill(v0->Pt());
1491 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_K0_Bin%d",ptbin)))->Fill(origV0);
a2ad7da1 1492 }
1493 N_Kaons++;
1494 } // end kaon case
1495
1496 //leading track correlations fill
1497 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1498 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0); //c and b D0
1499 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0); //c or b D0
1500 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);
1501 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);
1502 }
1503 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1504 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
1505 ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar); //c or b D0
1506 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);
1507 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);
1508 }
1509
1510 } //end fReadMC = 1
1511
1512 if (!fAlreadyFilled) {
1513 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Charg_Bin%d",ptbin)))->Fill(N_Charg);
1514 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Kcharg_Bin%d",ptbin)))->Fill(N_Kcharg);
f80e7bba 1515 ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_K0_Bin%d",ptbin)))->Fill(N_Kaons);
a2ad7da1 1516 }
1517
1518 fAlreadyFilled=kTRUE; //distribution plots for tracks filled
1519
1520}
1521
1522//_________________________________________________________________________________________________
1523Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1524 //
1525 // checks on particle (#) origin:
1526 // 0) Not HF
1527 // 1) D->#
1528 // 2) D->X->#
1529 // 3) B->#
1530 // 4) B->X-># (X!=D)
1531 // 5) B->D->#
1532 // 6) B->D->X->#
1533 // 7) c hadronization
1534 // 8) b hadronization
1535 //
1536 printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
1537
1538 Int_t pdgGranma = 0;
1539 Int_t mother = 0;
1540 mother = mcPartCandidate->GetMother();
1541 Int_t istep = 0;
1542 Int_t abspdgGranma =0;
1543 Bool_t isFromB=kFALSE;
1544 Bool_t isDdaugh=kFALSE;
1545 Bool_t isDchaindaugh=kFALSE;
1546 Bool_t isBdaugh=kFALSE;
1547 Bool_t isBchaindaugh=kFALSE;
1548 Bool_t isQuarkFound=kFALSE;
1549
1550 while (mother > 0){
1551 istep++;
1552 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1553 if (mcGranma){
1554 pdgGranma = mcGranma->GetPdgCode();
1555 abspdgGranma = TMath::Abs(pdgGranma);
1556 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1557 isBchaindaugh=kTRUE;
1558 if(istep==1) isBdaugh=kTRUE;
1559 }
1560 if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
1561 isDchaindaugh=kTRUE;
1562 if(istep==1) isDdaugh=kTRUE;
1563 }
1564 if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
1565 mother = mcGranma->GetMother();
1566 }else{
1567 AliError("Failed casting the mother particle!");
1568 break;
1569 }
1570 }
1571
1572 //decides what to return based on the flag status
1573 if(isQuarkFound) {
1574 if(!isFromB) { //charm
1575 if(isDdaugh) return 1; //charm immediate
1576 else if(isDchaindaugh) return 2; //charm chain
1577 else return 7; //charm hadronization
1578 }
1579 else { //beauty
1580 if(isBdaugh) return 3; //b immediate
1581 else if(isBchaindaugh) { //b chain
1582 if(isDchaindaugh) {
1583 if(isDdaugh) return 5; //d immediate
1584 return 6; //d chain
1585 }
1586 else return 4; //b, not d
1587 }
1588 else return 8; //b hadronization
1589 }
1590 }
1591 else return 0; //no HF quark
1592}
1593
1594
1595//________________________________________________________________________
1596Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
1597 //
1598 //give the pt bin where the pt lies.
1599 //
1600 Int_t ptbin=-1;
1601 if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
1602
1603 Int_t i = 0;
1604 while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
1605
1606 return ptbin;
1607}
1608
1609//________________________________________________________________________
1610Double_t AliAnalysisTaskSED0Correlations::PtWeig(Int_t ptbin, Double_t x) const {
1611 //
1612 //gives the weight for Weighted Corrs (inverse of pT distribution, from 3 D0 pt bins)
1613 //
1614 if(x>11.5) x=11.5; if(x<0.3) x=0.3; //bounds for fit of distributions!
1615 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)));
1616 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)));
1617 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)));
1618
1619 return 0; //for other bins plot is disabled!
1620}
1621
1622//________________________________________________________________________
1623void AliAnalysisTaskSED0Correlations::GetImpParameter(AliAODTrack *track, AliAODEvent* aod, Double_t &d0, Double_t &d0err) const {
1624 //
1625 //calculates d0 and error on d0 for the track
1626 //
1627 Double_t d0z0[2],covd0z0[3];
1628 AliESDtrack esdTrack(track); // AliESDTrack conversion of AOD track
1629 esdTrack.PropagateToDCA(aod->GetPrimaryVertex(),aod->GetMagneticField(),10000.,d0z0,covd0z0);
1630 d0 = TMath::Abs(d0z0[0]); // impact parameter
1631 d0err = TMath::Sqrt(covd0z0[0]); // resolution of impact parameter
1632}
1633
1634//________________________________________________________________________
1635Bool_t AliAnalysisTaskSED0Correlations::TrackSelectedInLoop(AliAODTrack* track, AliAODRecoDecayHF2Prong *d, AliAODEvent *aod, Int_t ptbin, Int_t idDaughs[]) const {
1636 //
1637 //rejection of tracks in loop for correlations
1638 //
1639 Bool_t output = kTRUE;
1640 AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
1641 Double_t bz = aod->GetMagneticField();
f80e7bba 1642 Double_t mD0, mD0bar;
1643 d->InvMassD0(mD0,mD0bar);
1644
a2ad7da1 1645 if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1] || track->GetID() < 0) output = kFALSE; //discards daughters of candidate
1646 if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) output = kFALSE; //discard tracks outside pt range for hadrons/K
1647 if(!fCutsTracks->IsHadronSelected(track,vtx,bz)) { //track discarded by quality cuts
1648 ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(1);
1649 output = kFALSE;
1650 } else ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
1651 if(!fCutsTracks->InvMassDstarRejection(d,track,fIsSelectedCandidate)) {
f80e7bba 1652 if (fReadMC == 0 && (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3)) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1653 if (fReadMC == 0 && fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1654 if (fReadMC == 1) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,1.864);
a2ad7da1 1655 output = kFALSE;
f80e7bba 1656 } else {
1657 if (fReadMC == 0 && (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3)) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1658 if (fReadMC == 0 && fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1659 if (fReadMC == 1) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,1.864);
1660 }
a2ad7da1 1661
1662 return output;
1663}
1664
1665//---------------------------------------------------------------------------
1666Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
1667{
1668 //
1669 // Selection for K0 hypotheses
1670 // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
1671 // 2 = no previous selections
1672
1673 if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
1674
1675 AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
1676 AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
1677
1678 if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
1679 if(v0->Pt() < 3. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0040) return kFALSE;
1680 if(v0->Pt() > 3. && v0->Pt() < 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0052) return kFALSE;
1681 if(v0->Pt() > 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0075) return kFALSE;
1682 }
1683
1684 //This part removes double counting for swapped tracks!
1685 Int_t i = 0; //while loop (until the last-written entry pair of ID!
1686 while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
1687 if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
1688 (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
1689 i++;
1690 }
1691 idArrayV0[i][0]=v0Daug1->GetID();
1692 idArrayV0[i][1]=v0Daug2->GetID();
1693
1694 return kTRUE;
1695}
1696
1697//________________________________________________________________________
1698void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
1699
1700 cout << "--------------------------\n";
1701 cout << "PtBins = " << fNPtBinsCorr << "\n";
1702 cout << "PtBin limits--------------\n";
1703 for (int i=0; i<fNPtBinsCorr; i++) {
1704 cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
1705 }
1706 cout << "\n--------------------------\n";
1707 cout << "PtBin tresh. tracks low---\n";
1708 for (int i=0; i<fNPtBinsCorr; i++) {
1709 cout << fPtThreshLow.at(i) << "; ";
1710 }
1711 cout << "PtBin tresh. tracks up----\n";
1712 for (int i=0; i<fNPtBinsCorr; i++) {
1713 cout << fPtThreshUp.at(i) << "; ";
1714 }
1715 cout << "\n--------------------------\n";
1716 cout << "MC Truth = "<<fReadMC<<"\n";
1717 cout << "--------------------------\n";
1718}
1719