1 /**************************************************************************
2 * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // AliAnalysisTaskSE for D0 candidates (2Prongs)
21 // and hadrons correlations
23 // Authors: Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
24 // Fabio Colamaria, fabio.colamaria@ba.infn.it (correlations)
25 /////////////////////////////////////////////////////////////
27 #include <Riostream.h>
28 #include <TClonesArray.h>
35 #include <THnSparse.h>
36 #include <TDatabasePDG.h>
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"
59 ClassImp(AliAnalysisTaskSED0Correlations)
62 //________________________________________________________________________
63 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations():
71 fAlreadyFilled(kFALSE),
89 fIsSelectedCandidate(0),
92 fIsRejectSDDClusters(0),
95 // Default constructor
99 //________________________________________________________________________
100 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *name,AliRDHFCutsD0toKpi* cutsD0, AliHFAssociatedTrackCuts* cutsTrk):
101 AliAnalysisTaskSE(name),
108 fAlreadyFilled(kFALSE),
114 fCutsTracks(cutsTrk),
126 fIsSelectedCandidate(0),
129 fIsRejectSDDClusters(0),
132 // Default constructor
134 fNPtBins=cutsD0->GetNPtBins();
138 // Output slot #1 writes into a TList container (mass with cuts)
139 DefineOutput(1,TList::Class()); //My private output
140 // Output slot #2 writes into a TH1F container (number of events)
141 DefineOutput(2,TH1F::Class()); //My private output
142 // Output slot #3 writes into a AliRDHFD0toKpi container (cuts)
143 DefineOutput(3,AliRDHFCutsD0toKpi::Class()); //My private output
144 // Output slot #4 writes Normalization Counter
145 DefineOutput(4,AliNormalizationCounter::Class());
146 // Output slot #5 writes into a TList container (correl output)
147 DefineOutput(5,TList::Class()); //My private output
148 // Output slot #6 writes into a TList container (correl advanced)
149 DefineOutput(6,TList::Class()); //My private output
150 // Output slot #7 writes into a AliHFAssociatedTrackCuts container (cuts)
151 DefineOutput(7,AliHFAssociatedTrackCuts::Class()); //My private output
154 //________________________________________________________________________
155 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalysisTaskSED0Correlations &source):
156 AliAnalysisTaskSE(source),
157 fNPtBinsCorr(source.fNPtBinsCorr),
158 fBinLimsCorr(source.fBinLimsCorr),
159 fPtThreshLow(source.fPtThreshLow),
160 fPtThreshUp(source.fPtThreshUp),
161 fD0Eff(source.fD0Eff),
162 fEvents(source.fEvents),
163 fAlreadyFilled(source.fAlreadyFilled),
164 fOutputMass(source.fOutputMass),
165 fOutputCorr(source.fOutputCorr),
166 fOutputStudy(source.fOutputStudy),
167 fNentries(source.fNentries),
168 fCutsD0(source.fCutsD0),
169 fCutsTracks(source.fCutsTracks),
170 fCorrelatorTr(source.fCorrelatorTr),
171 fCorrelatorKc(source.fCorrelatorKc),
172 fCorrelatorK0(source.fCorrelatorK0),
173 fReadMC(source.fReadMC),
174 fRecoTr(source.fRecoTr),
175 fRecoD0(source.fRecoD0),
176 fSelEvType(source.fSelEvType),
177 fMixing(source.fMixing),
178 fCounter(source.fCounter),
179 fNPtBins(source.fNPtBins),
180 fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
181 fIsSelectedCandidate(source.fIsSelectedCandidate),
183 fEtaForCorrel(source.fEtaForCorrel),
184 fIsRejectSDDClusters(source.fIsRejectSDDClusters),
185 fFillGlobal(source.fFillGlobal)
190 //________________________________________________________________________
191 AliAnalysisTaskSED0Correlations::~AliAnalysisTaskSED0Correlations()
214 delete fCorrelatorTr;
218 delete fCorrelatorKc;
222 delete fCorrelatorK0;
231 //______________________________________________________________________________
232 AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(const AliAnalysisTaskSED0Correlations& orig)
235 if (&orig == this) return *this; //if address is the same (same object), returns itself
237 AliAnalysisTaskSE::operator=(orig); //Uses the AliAnalysisTaskSE operator to assign the inherited part of the class
238 fNPtBinsCorr = orig.fNPtBinsCorr;
239 fBinLimsCorr = orig.fBinLimsCorr;
240 fPtThreshLow = orig.fPtThreshLow;
241 fPtThreshUp = orig.fPtThreshUp;
242 fD0Eff = orig.fD0Eff;
243 fEvents = orig.fEvents;
244 fAlreadyFilled = orig.fAlreadyFilled;
245 fOutputMass = orig.fOutputMass;
246 fOutputCorr = orig.fOutputCorr;
247 fOutputStudy = orig.fOutputStudy;
248 fNentries = orig.fNentries;
249 fCutsD0 = orig.fCutsD0;
250 fCutsTracks = orig.fCutsTracks;
251 fCorrelatorTr = orig.fCorrelatorTr;
252 fCorrelatorKc = orig.fCorrelatorKc;
253 fCorrelatorK0 = orig.fCorrelatorK0;
254 fReadMC = orig.fReadMC;
255 fRecoTr = orig.fRecoTr;
256 fRecoD0 = orig.fRecoD0;
257 fSelEvType = orig.fSelEvType;
258 fMixing = orig.fMixing;
259 fCounter = orig.fCounter;
260 fNPtBins = orig.fNPtBins;
261 fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
262 fIsSelectedCandidate = orig.fIsSelectedCandidate;
264 fEtaForCorrel = orig.fEtaForCorrel;
265 fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
266 fFillGlobal = orig.fFillGlobal;
268 return *this; //returns pointer of the class
271 //________________________________________________________________________
272 void AliAnalysisTaskSED0Correlations::Init()
276 if(fDebug > 1) printf("AnalysisTaskSED0Correlations::Init() \n");
278 //Copy of cuts objects
279 AliRDHFCutsD0toKpi* copyfCutsD0 = new AliRDHFCutsD0toKpi(*fCutsD0);
280 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
281 copyfCutsD0->SetName(nameoutput);
283 //needer to clear completely the objects inside with Clear() method
285 PostData(3,copyfCutsD0);
286 PostData(7,fCutsTracks);
291 //________________________________________________________________________
292 void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
295 // Create the output container
297 if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
299 //HFCorrelator creation and definition
300 fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys);
301 fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys);
302 fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys);
303 fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);// set the Delta Phi Interval you want
304 fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);
305 fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2+TMath::Pi()/32.,3*TMath::Pi()/2+TMath::Pi()/32.);
306 fCorrelatorTr->SetEventMixing(fMixing);// sets the analysis on a single event (kFALSE) or mixed events (kTRUE)
307 fCorrelatorKc->SetEventMixing(fMixing);
308 fCorrelatorK0->SetEventMixing(fMixing);
309 fCorrelatorTr->SetAssociatedParticleType(1);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
310 fCorrelatorKc->SetAssociatedParticleType(2);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
311 fCorrelatorK0->SetAssociatedParticleType(3);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
312 fCorrelatorTr->SetApplyDisplacementCut(2); //0: don't calculate d0; 1: return d0; 2: return d0/d0err
313 fCorrelatorKc->SetApplyDisplacementCut(2);
314 fCorrelatorK0->SetApplyDisplacementCut(0);
315 fCorrelatorTr->SetUseMC(fReadMC);// sets Montecarlo flag
316 fCorrelatorKc->SetUseMC(fReadMC);
317 fCorrelatorK0->SetUseMC(fReadMC);
318 fCorrelatorTr->SetUseReco(fRecoTr);// sets (if MC analysis) wheter to analyze Reco or Kinem tracks
319 fCorrelatorKc->SetUseReco(fRecoTr);
320 fCorrelatorK0->SetUseReco(fRecoTr);
321 fCorrelatorKc->SetPIDmode(2); //switch for K+/- PID option
322 Bool_t pooldefTr = fCorrelatorTr->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
323 Bool_t pooldefKc = fCorrelatorKc->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
324 Bool_t pooldefK0 = fCorrelatorK0->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
325 if(!pooldefTr) AliInfo("Warning:: Event pool not defined properly");
326 if(!pooldefKc) AliInfo("Warning:: Event pool not defined properly");
327 if(!pooldefK0) AliInfo("Warning:: Event pool not defined properly");
329 // Several histograms are more conveniently managed in a TList
330 fOutputMass = new TList();
331 fOutputMass->SetOwner();
332 fOutputMass->SetName("listMass");
334 fOutputCorr = new TList();
335 fOutputCorr->SetOwner();
336 fOutputCorr->SetName("correlationslist");
338 fOutputStudy = new TList();
339 fOutputStudy->SetOwner();
340 fOutputStudy->SetName("MCstudyplots");
342 TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassWg=" ",nameSgnWg=" ", nameBkgWg=" ", nameRflWg=" ";
344 for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
346 nameMass="histMass_";
348 nameMassWg="histMass_WeigD0Eff_";
352 nameSgnWg="histSgn_WeigD0Eff_";
356 nameBkgWg="histBkg_WeigD0Eff_";
360 nameRflWg="histRfl_WeigD0Eff_";
363 //histograms of invariant mass distributions
367 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
368 TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
372 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
373 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
374 TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
375 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
376 TH1F* tmpBtWg = new TH1F(nameBkgWg.Data(), "Background invariant mass - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
381 fOutputMass->Add(tmpSt);
382 fOutputMass->Add(tmpStWg);
383 fOutputMass->Add(tmpRt);
384 fOutputMass->Add(tmpRtWg);
385 fOutputMass->Add(tmpBt);
386 fOutputMass->Add(tmpBtWg);
390 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",120,1.5648,2.1648);
392 fOutputMass->Add(tmpMt);
393 //mass weighted by 1/D0eff
394 TH1F* tmpMtwg = new TH1F(nameMassWg.Data(),"D^{0} invariant mass - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
396 fOutputMass->Add(tmpMtwg);
399 const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
401 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", 20,-0.5,19.5);
403 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
404 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
405 fReadMC ? fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected") : fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
406 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
407 fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
408 fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
409 if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
410 if(fReadMC && fSys==0){
411 fNentries->GetXaxis()->SetBinLabel(12,"K");
412 fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
414 fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
415 fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
416 if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
417 if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
418 fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
419 fNentries->GetXaxis()->SetBinLabel(19,"nEventsSelected");
420 if(fReadMC) fNentries->GetXaxis()->SetBinLabel(20,"nEvsWithProdMech");
421 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
423 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
426 CreateCorrelationsObjs(); //creates histos for correlations analysis
429 PostData(1,fOutputMass);
430 PostData(2,fNentries);
431 PostData(4,fCounter);
432 PostData(5,fOutputCorr);
433 PostData(6,fOutputStudy);
438 //________________________________________________________________________
439 void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
441 // Execute analysis for current event:
442 // heavy flavor candidates association to MC truth
443 //cout<<"I'm in UserExec"<<endl;
447 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
448 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
449 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
450 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
451 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
452 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
453 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
454 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
455 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
457 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
460 TString bname="D0toKpi";
462 TClonesArray *inputArray=0;
464 if(!aod && AODEvent() && IsStandardAOD()) {
465 // In case there is an AOD handler writing a standard AOD, use the AOD
466 // event in memory rather than the input (ESD) event.
467 aod = dynamic_cast<AliAODEvent*> (AODEvent());
468 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
469 // have to taken from the AOD event hold by the AliAODExtension
470 AliAODHandler* aodHandler = (AliAODHandler*)
471 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
473 if(aodHandler->GetExtensions()) {
474 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
475 AliAODEvent* aodFromExt = ext->GetAOD();
476 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
479 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
482 if(!inputArray || !aod) {
483 printf("AliAnalysisTaskSED0Correlations::UserExec: input branch not found!\n");
487 // fix for temporary bug in ESDfilter
488 // the AODs with null vertex pointer didn't pass the PhysSel
489 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
491 TClonesArray *mcArray = 0;
492 AliAODMCHeader *mcHeader = 0;
496 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
498 printf("AliAnalysisTaskSED0Correlations::UserExec: MC particles branch not found!\n");
503 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
505 printf("AliAnalysisTaskSED0Correlations::UserExec: MC header branch not found!\n");
510 //histogram filled with 1 for every AOD
512 fCounter->StoreEvent(aod,fCutsD0,fReadMC);
514 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
515 TString trigclass=aod->GetFiredTriggerClasses();
516 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
518 if(!fCutsD0->IsEventSelected(aod)) {
519 if(fCutsD0->GetWhyRejection()==1) // rejected for pileup
521 if(fSys==1 && (fCutsD0->GetWhyRejection()==2 || fCutsD0->GetWhyRejection()==3)) fNentries->Fill(15);
522 if(fCutsD0->GetWhyRejection()==7) fNentries->Fill(17);
526 fNentries->Fill(18); //event selected after selection
528 //Setting PIDResponse for associated tracks
529 fCorrelatorTr->SetPidAssociated();
530 fCorrelatorKc->SetPidAssociated();
531 fCorrelatorK0->SetPidAssociated();
533 //Selection on production type (MC)
534 if(fReadMC && fSelEvType){
536 Bool_t isMCeventgood = kFALSE;
538 Int_t eventType = mcHeader->GetEventType();
539 Int_t NMCevents = fCutsTracks->GetNofMCEventType();
541 for(Int_t k=0; k<NMCevents; k++){
542 Int_t * MCEventType = fCutsTracks->GetMCEventType();
544 if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
545 ((TH1D*)fOutputStudy->FindObject("EventTypeMC"))->Fill(eventType);
548 if(NMCevents && !isMCeventgood){
549 std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
552 fNentries->Fill(19); //event with particular production type
557 // Check the Nb of SDD clusters
558 if (fIsRejectSDDClusters) {
559 Bool_t skipEvent = kFALSE;
561 if (aod) ntracks = aod->GetNTracks();
562 for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
564 AliAODTrack * track = aod->GetTrack(itrack);
565 if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
571 if (skipEvent) return;
574 //HFCorrelators initialization (for this event)
575 fCorrelatorTr->SetAODEvent(aod); // set the AOD event from which you are processing
576 fCorrelatorKc->SetAODEvent(aod);
577 fCorrelatorK0->SetAODEvent(aod);
578 Bool_t correlatorONTr = fCorrelatorTr->Initialize(); // initialize the pool for event mixing
579 Bool_t correlatorONKc = fCorrelatorKc->Initialize();
580 Bool_t correlatorONK0 = fCorrelatorK0->Initialize();
581 if(!correlatorONTr) {AliInfo("AliHFCorrelator (tracks) didn't initialize the pool correctly or processed a bad event"); return;}
582 if(!correlatorONKc) {AliInfo("AliHFCorrelator (charged K) didn't initialize the pool correctly or processed a bad event"); return;}
583 if(!correlatorONK0) {AliInfo("AliHFCorrelator (K0) didn't initialize the pool correctly or processed a bad event"); return;}
586 fCorrelatorTr->SetMCArray(mcArray); // set the TClonesArray *fmcArray for analysis on monte carlo
587 fCorrelatorKc->SetMCArray(mcArray);
588 fCorrelatorK0->SetMCArray(mcArray);
591 // AOD primary vertex
592 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
595 TString primTitle = vtx1->GetTitle();
596 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
600 //Reset flag for tracks distributions fill
601 fAlreadyFilled=kFALSE;
603 //***** Loop over D0 candidates *****
604 Int_t nInD0toKpi = inputArray->GetEntriesFast();
605 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
607 if(fFillGlobal) { //loop on V0 and tracks for each event, to fill Pt distr. and InvMass distr.
609 TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
610 Int_t pdgCodes[2] = {211,211};
611 Int_t idArrayV0[v0array->GetEntriesFast()][2];
612 for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
613 for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
614 AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
615 if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
616 if(fReadMC && fRecoTr && (v0->MatchToMC(310,mcArray,2,pdgCodes)<0)) continue; //310 = K0s, 311 = K0 generico!!
617 ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
618 ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
622 for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tacks
623 AliAODTrack * track = aod->GetTrack(itrack);
624 //rejection of tracks
625 if(track->GetID() < 0) continue; //discard negative ID tracks
626 if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
627 if(!fCutsTracks->IsHadronSelected(track) || !fCutsTracks->CheckHadronKinematic(track->Pt(),0.1)) continue; //0.1 = dummy (d0 value, no cut on it for me)
628 //pT distribution (in all events), charg and hadr cases
629 ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt());
630 if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
633 } //end of loops for global plot fill
635 Int_t nSelectedloose=0,nSelectedtight=0;
637 //RecoD0 case ************************************************
640 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
641 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
643 if(d->Pt()<2.) continue; //to save time and merging memory...
645 if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
647 continue; //skip the D0 from Dstar
650 if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
654 if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
656 Int_t ptbin=fCutsD0->PtBin(d->Pt());
657 if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
659 fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
660 if(!fIsSelectedCandidate) continue;
663 Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(d->Phi());
664 phiD0 = fCorrelatorKc->SetCorrectPhiRange(d->Phi()); //bad usage, but returns a Double_t...
665 phiD0 = fCorrelatorK0->SetCorrectPhiRange(d->Phi());
666 fCorrelatorTr->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta()); // sets the parameters of the trigger particles that are needed
667 fCorrelatorKc->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
668 fCorrelatorK0->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
669 fCorrelatorTr->SetD0Properties(d,fIsSelectedCandidate); //sets special properties for D0
670 fCorrelatorKc->SetD0Properties(d,fIsSelectedCandidate);
671 fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
674 if (TMath::Abs(d->Eta())<fEtaForCorrel) CalculateCorrelations(d); //correlations on real data
675 } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
676 if (TMath::Abs(d->Eta())<fEtaForCorrel) {
677 Int_t pdgDgD0toKpi[2]={321,211};
678 Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
679 if (labD0>-1) CalculateCorrelations(d,labD0,mcArray);
683 FillMassHists(d,mcArray,fCutsD0,fOutputMass);
687 //End RecoD0 case ************************************************
689 //MCKineD0 case ************************************************
690 if(fReadMC && !fRecoD0) {
692 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { //Loop over all the tracks of MCArray
693 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
695 AliWarning("Particle not found in tree, skipping");
699 if(TMath::Abs(mcPart->GetPdgCode()) == 421){ // THIS IS A D0
700 if (fCutsD0->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y()) ) {
703 if(fSys==0) fNentries->Fill(6);
704 Int_t ptbin=fCutsD0->PtBin(mcPart->Pt());
705 if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
708 Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(mcPart->Phi());
709 phiD0 = fCorrelatorKc->SetCorrectPhiRange(mcPart->Phi()); //bad usage, but returns a Double_t...
710 phiD0 = fCorrelatorK0->SetCorrectPhiRange(mcPart->Phi());
711 fCorrelatorTr->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta()); // sets the parameters of the trigger particles that are needed
712 fCorrelatorKc->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
713 fCorrelatorK0->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
714 //fCorrelatorTr->SetD0Properties(mcPart,fIsSelectedCandidate); //needed for D* soft pions rejection, useless in MCKine
715 //fCorrelatorKc->SetD0Properties(mcPart,fIsSelectedCandidate);
716 //fCorrelatorK0->SetD0Properties(mcPart,fIsSelectedCandidate);
718 if (TMath::Abs(mcPart->Eta())<fEtaForCorrel) {
720 //Removal of D0 from D* feeddown! This solves also the problem of soft pions, now excluded
721 /* Int_t mother = mcPart->GetMother();
722 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
723 if(!mcMoth) continue;
724 if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;
726 if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
727 else fIsSelectedCandidate = 2;
729 TString fillthis="histSgn_"; fillthis+=ptbin;
730 ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
732 CalculateCorrelationsMCKine(mcPart,mcArray);
738 //End MCKineD0 case ************************************************
740 if(fMixing /* && fAlreadyFilled*/) { // update the pool for Event Mixing, if: enabled, event is ok, at least a SelD0 found! (fAlreadyFilled's role!)
741 Bool_t updatedTr = fCorrelatorTr->PoolUpdate();
742 Bool_t updatedKc = fCorrelatorKc->PoolUpdate();
743 Bool_t updatedK0 = fCorrelatorK0->PoolUpdate();
744 if(!updatedTr || !updatedKc || !updatedK0) AliInfo("Pool was not updated");
747 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
748 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
751 PostData(1,fOutputMass);
752 PostData(2,fNentries);
753 PostData(4,fCounter);
754 PostData(5,fOutputCorr);
755 PostData(6,fOutputStudy);
760 //____________________________________________________________________________
761 void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
763 // function used in UserExec to fill mass histograms:
765 if (!fIsSelectedCandidate) return;
767 if(fDebug>2) cout<<"Candidate selected"<<endl;
769 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
770 Int_t ptbin = cuts->PtBin(part->Pt());
773 Int_t pdgDgD0toKpi[2]={321,211};
775 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)
777 //count candidates selected by cuts
779 //count true D0 selected by cuts
780 if (fReadMC && labD0>=0) fNentries->Fill(2);
782 if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
786 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
787 Int_t pdgD0 = partD0->GetPdgCode();
788 if (pdgD0==421){ //D0
791 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
792 fillthis="histSgn_WeigD0Eff_";
794 Double_t effD0 = fD0Eff.at(ptbin);
795 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
796 } else{ //it was a D0bar
799 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
800 fillthis="histRfl_WeigD0Eff_";
802 Double_t effD0 = fD0Eff.at(ptbin);
803 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
808 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
809 fillthis="histBkg_WeigD0Eff_";
811 Double_t effD0 = fD0Eff.at(ptbin);
812 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
815 fillthis="histMass_";
817 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
818 fillthis="histMass_WeigD0Eff_";
820 Double_t effD0 = fD0Eff.at(ptbin);
821 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
825 if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
829 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
830 Int_t pdgD0 = partD0->GetPdgCode();
832 if (pdgD0==-421){ //D0bar
835 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
836 fillthis="histSgn_WeigD0Eff_";
838 Double_t effD0 = fD0Eff.at(ptbin);
839 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
843 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
844 fillthis="histRfl_WeigD0Eff_";
846 Double_t effD0 = fD0Eff.at(ptbin);
847 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
849 } else {//background or LS
852 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
853 fillthis="histBkg_WeigD0Eff_";
855 Double_t effD0 = fD0Eff.at(ptbin);
856 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
859 fillthis="histMass_";
861 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
862 fillthis="histMass_WeigD0Eff_";
864 Double_t effD0 = fD0Eff.at(ptbin);
865 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
872 //________________________________________________________________________
873 void AliAnalysisTaskSED0Correlations::Terminate(Option_t */*option*/)
875 // Terminate analysis
877 if(fDebug > 1) printf("AnalysisTaskSED0Correlations: Terminate() \n");
879 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
881 printf("ERROR: fOutputMass not available\n");
885 fNentries = dynamic_cast<TH1F*>(GetOutputData(2));
888 printf("ERROR: fNEntries not available\n");
892 fCutsD0 = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(3));
894 printf("ERROR: fCuts not available\n");
898 fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(4));
900 printf("ERROR: fCounter not available\n");
903 fOutputCorr = dynamic_cast<TList*> (GetOutputData(5));
905 printf("ERROR: fOutputCorr not available\n");
908 fOutputStudy = dynamic_cast<TList*> (GetOutputData(6));
910 printf("ERROR: fOutputStudy not available\n");
913 fCutsTracks = dynamic_cast<AliHFAssociatedTrackCuts*>(GetOutputData(7));
915 printf("ERROR: fCutsTracks not available\n");
922 //_________________________________________________________________________________________________
923 Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
925 // checking whether the mother of the particles come from a charm or a bottom quark
927 printf("AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
931 mother = mcPartCandidate->GetMother();
932 Int_t abspdgGranma =0;
933 Bool_t isFromB=kFALSE;
934 Bool_t isQuarkFound=kFALSE;
937 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
939 pdgGranma = mcMoth->GetPdgCode();
940 abspdgGranma = TMath::Abs(pdgGranma);
941 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
944 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
945 mother = mcMoth->GetMother();
947 AliError("Failed casting the mother particle!");
953 if(isFromB) return 5;
959 //________________________________________________________________________
960 void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
963 TString namePlot = "";
965 //These for limits in THnSparse (one per bin, same limits).
966 //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
967 Int_t nBinsPhi[5] = {32,150,6,3,16};
968 Double_t binMinPhi[5] = {-TMath::Pi()/2.+TMath::Pi()/32.,1.6,0.,0.,-1.6}; //is the minimum for all the bins
969 Double_t binMaxPhi[5] = {3*TMath::Pi()/2.+TMath::Pi()/32.,2.2,3.0,3.,1.6}; //is the maximum for all the bins
971 //Vars: DeltaPhi, InvMass, DeltaEta
972 Int_t nBinsMix[3] = {32,150,16};
973 Double_t binMinMix[3] = {-TMath::Pi()/2.+TMath::Pi()/32.,1.6,-1.6}; //is the minimum for all the bins
974 Double_t binMaxMix[3] = {3*TMath::Pi()/2.+TMath::Pi()/32.,2.2,1.6}; //is the maximum for all the bins
976 for(Int_t i=0;i<fNPtBinsCorr;i++){
979 //THnSparse plots: correlations for various invariant mass (MC and data)
980 namePlot="hPhi_K0_Bin";
983 THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
985 fOutputCorr->Add(hPhiK);
987 namePlot="hPhi_Kcharg_Bin";
990 THnSparseF *hPhiH = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
992 fOutputCorr->Add(hPhiH);
994 namePlot="hPhi_Charg_Bin";
997 THnSparseF *hPhiC = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
999 fOutputCorr->Add(hPhiC);
1001 //histos for c/b origin for D0 (MC only)
1004 //generic origin for tracks
1005 namePlot="hPhi_K0_From_c_Bin";
1008 THnSparseF *hPhiK_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1010 fOutputCorr->Add(hPhiK_c);
1012 namePlot="hPhi_Kcharg_From_c_Bin";
1015 THnSparseF *hPhiH_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1017 fOutputCorr->Add(hPhiH_c);
1019 namePlot="hPhi_Charg_From_c_Bin";
1022 THnSparseF *hPhiC_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1024 fOutputCorr->Add(hPhiC_c);
1026 namePlot="hPhi_K0_From_b_Bin";
1029 THnSparseF *hPhiK_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1031 fOutputCorr->Add(hPhiK_b);
1033 namePlot="hPhi_Kcharg_From_b_Bin";
1036 THnSparseF *hPhiH_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1038 fOutputCorr->Add(hPhiH_b);
1040 namePlot="hPhi_Charg_From_b_Bin";
1043 THnSparseF *hPhiC_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1045 fOutputCorr->Add(hPhiC_b);
1047 //HF-only tracks (c for c->D0, b for b->D0)
1048 namePlot="hPhi_K0_HF_From_c_Bin";
1051 THnSparseF *hPhiK_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1052 hPhiK_HF_c->Sumw2();
1053 fOutputCorr->Add(hPhiK_HF_c);
1055 namePlot="hPhi_Kcharg_HF_From_c_Bin";
1058 THnSparseF *hPhiH_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1059 hPhiH_HF_c->Sumw2();
1060 fOutputCorr->Add(hPhiH_HF_c);
1062 namePlot="hPhi_Charg_HF_From_c_Bin";
1065 THnSparseF *hPhiC_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1066 hPhiC_HF_c->Sumw2();
1067 fOutputCorr->Add(hPhiC_HF_c);
1069 namePlot="hPhi_K0_HF_From_b_Bin";
1072 THnSparseF *hPhiK_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1073 hPhiK_HF_b->Sumw2();
1074 fOutputCorr->Add(hPhiK_HF_b);
1076 namePlot="hPhi_Kcharg_HF_From_b_Bin";
1079 THnSparseF *hPhiH_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1080 hPhiH_HF_b->Sumw2();
1081 fOutputCorr->Add(hPhiH_HF_b);
1083 namePlot="hPhi_Charg_HF_From_b_Bin";
1086 THnSparseF *hPhiC_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1087 hPhiC_HF_b->Sumw2();
1088 fOutputCorr->Add(hPhiC_HF_b);
1090 namePlot="hPhi_K0_NonHF_Bin";
1093 THnSparseF *hPhiK_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1095 fOutputCorr->Add(hPhiK_Non);
1097 namePlot="hPhi_Kcharg_NonHF_Bin";
1100 THnSparseF *hPhiH_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1102 fOutputCorr->Add(hPhiH_Non);
1104 namePlot="hPhi_Charg_NonHF_Bin";
1107 THnSparseF *hPhiC_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1109 fOutputCorr->Add(hPhiC_Non);
1112 //leading hadron correlations
1113 namePlot="hPhi_Lead_Bin";
1116 THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1118 fOutputCorr->Add(hCorrLead);
1121 namePlot="hPhi_Lead_From_c_Bin";
1124 THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1125 hCorrLead_c->Sumw2();
1126 fOutputCorr->Add(hCorrLead_c);
1128 namePlot="hPhi_Lead_From_b_Bin";
1131 THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1132 hCorrLead_b->Sumw2();
1133 fOutputCorr->Add(hCorrLead_b);
1135 namePlot="hPhi_Lead_HF_From_c_Bin";
1138 THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1139 hCorrLead_HF_c->Sumw2();
1140 fOutputCorr->Add(hCorrLead_HF_c);
1142 namePlot="hPhi_Lead_HF_From_b_Bin";
1145 THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1146 hCorrLead_HF_b->Sumw2();
1147 fOutputCorr->Add(hCorrLead_HF_b);
1149 namePlot="hPhi_Lead_NonHF_Bin";
1152 THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1153 hCorrLead_Non->Sumw2();
1154 fOutputCorr->Add(hCorrLead_Non);
1157 //pT weighted correlations
1158 namePlot="hPhi_Weig_Bin";
1161 THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1162 fOutputCorr->Add(hCorrWeig);
1165 namePlot="hPhi_Weig_From_c_Bin";
1168 THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1169 fOutputCorr->Add(hCorrWeig_c);
1171 namePlot="hPhi_Weig_From_b_Bin";
1174 THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1175 fOutputCorr->Add(hCorrWeig_b);
1177 namePlot="hPhi_Weig_HF_From_c_Bin";
1180 THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1181 fOutputCorr->Add(hCorrWeig_HF_c);
1183 namePlot="hPhi_Weig_HF_From_b_Bin";
1186 THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1187 fOutputCorr->Add(hCorrWeig_HF_b);
1189 namePlot="hPhi_Weig_NonHF_Bin";
1192 THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1193 fOutputCorr->Add(hCorrWeig_Non);
1196 //pT distribution histos
1197 namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
1198 TH1F *hPtC = new TH1F(namePlot.Data(), "Charged track pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1199 hPtC->SetMinimum(0);
1200 fOutputStudy->Add(hPtC);
1202 namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
1203 TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1204 hPtH->SetMinimum(0);
1205 fOutputStudy->Add(hPtH);
1207 namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
1208 TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1209 hPtK->SetMinimum(0);
1210 fOutputStudy->Add(hPtK);
1212 //D* feeddown pions rejection histos
1213 namePlot = "hDstarPions_Bin"; namePlot+=i;
1214 TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,150,1.6,2.2);
1215 hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
1216 hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
1217 hDstarPions->SetMinimum(0);
1218 fOutputStudy->Add(hDstarPions);
1223 //THnSparse plots for event mixing!
1224 namePlot="hPhi_K0_Bin";
1225 namePlot+=i;namePlot+="_EvMix";
1227 THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1228 hPhiK_EvMix->Sumw2();
1229 fOutputCorr->Add(hPhiK_EvMix);
1231 namePlot="hPhi_Kcharg_Bin";
1232 namePlot+=i;namePlot+="_EvMix";
1234 THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1235 hPhiH_EvMix->Sumw2();
1236 fOutputCorr->Add(hPhiH_EvMix);
1238 namePlot="hPhi_Charg_Bin";
1239 namePlot+=i;namePlot+="_EvMix";
1241 THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1242 hPhiC_EvMix->Sumw2();
1243 fOutputCorr->Add(hPhiC_EvMix);
1249 TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1250 hCountC->SetMinimum(0);
1251 fOutputStudy->Add(hCountC);
1253 TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1254 hCountH->SetMinimum(0);
1255 fOutputStudy->Add(hCountH);
1257 TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1258 hCountK->SetMinimum(0);
1259 fOutputStudy->Add(hCountK);
1263 TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1264 fOutputStudy->Add(hEventTypeMC);
1267 if (fFillGlobal) { //all-events plots
1269 TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1270 hPtCAll->SetMinimum(0);
1271 fOutputStudy->Add(hPtCAll);
1273 TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1274 hPtHAll->SetMinimum(0);
1275 fOutputStudy->Add(hPtHAll);
1277 TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1278 hPtKAll->SetMinimum(0);
1279 fOutputStudy->Add(hPtKAll);
1283 TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1284 hPhiDistCAll->SetMinimum(0);
1285 fOutputStudy->Add(hPhiDistCAll);
1287 TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1288 hPhiDistHAll->SetMinimum(0);
1289 fOutputStudy->Add(hPhiDistHAll);
1291 TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1292 hPhiDistKAll->SetMinimum(0);
1293 fOutputStudy->Add(hPhiDistKAll);
1295 TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1296 hPhiDistDAll->SetMinimum(0);
1297 fOutputStudy->Add(hPhiDistDAll);
1300 //K0 Invariant Mass plots
1301 TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1302 hK0MassInv->SetMinimum(0);
1303 fOutputStudy->Add(hK0MassInv);
1306 //for MC analysis only
1307 for(Int_t i=0;i<fNPtBinsCorr;i++) {
1309 if (fReadMC && !fMixing) {
1311 //displacement histos
1312 namePlot="histDispl_K0_Bin"; namePlot+=i;
1313 TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1314 hDisplK->SetMinimum(0);
1315 fOutputStudy->Add(hDisplK);
1317 namePlot="histDispl_K0_HF_Bin"; namePlot+=i;
1318 TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1319 hDisplK_HF->SetMinimum(0);
1320 fOutputStudy->Add(hDisplK_HF);
1322 namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1323 TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1324 hDisplHadr->SetMinimum(0);
1325 fOutputStudy->Add(hDisplHadr);
1327 namePlot="histDispl_Kcharg_HF_Bin"; namePlot+=i;
1328 TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1329 hDisplHadr_HF->SetMinimum(0);
1330 fOutputStudy->Add(hDisplHadr_HF);
1332 namePlot="histDispl_Charg_Bin"; namePlot+=i;
1333 TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1334 hDisplCharg->SetMinimum(0);
1335 fOutputStudy->Add(hDisplCharg);
1337 namePlot="histDispl_Charg_HF_Bin"; namePlot+=i;
1338 TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1339 hDisplCharg_HF->SetMinimum(0);
1340 fOutputStudy->Add(hDisplCharg_HF);
1342 namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
1343 TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1344 hDisplK_c->SetMinimum(0);
1345 fOutputStudy->Add(hDisplK_c);
1347 namePlot="histDispl_K0_HF_From_c_Bin"; namePlot+=i;
1348 TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1349 hDisplK_HF_c->SetMinimum(0);
1350 fOutputStudy->Add(hDisplK_HF_c);
1352 namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1353 TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1354 hDisplHadr_c->SetMinimum(0);
1355 fOutputStudy->Add(hDisplHadr_c);
1357 namePlot="histDispl_Kcharg_HF_From_c_Bin"; namePlot+=i;
1358 TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1359 hDisplHadr_HF_c->SetMinimum(0);
1360 fOutputStudy->Add(hDisplHadr_HF_c);
1362 namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1363 TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1364 hDisplCharg_c->Sumw2();
1365 hDisplCharg_c->SetMinimum(0);
1366 fOutputStudy->Add(hDisplCharg_c);
1368 namePlot="histDispl_Charg_HF_From_c_Bin"; namePlot+=i;
1369 TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1370 hDisplCharg_HF_c->SetMinimum(0);
1371 fOutputStudy->Add(hDisplCharg_HF_c);
1373 namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
1374 TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1375 hDisplK_b->SetMinimum(0);
1376 fOutputStudy->Add(hDisplK_b);
1378 namePlot="histDispl_K0_HF_From_b_Bin"; namePlot+=i;
1379 TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1380 hDisplK_HF_b->SetMinimum(0);
1381 fOutputStudy->Add(hDisplK_HF_b);
1383 namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1384 TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1385 hDisplHadr_b->SetMinimum(0);
1386 fOutputStudy->Add(hDisplHadr_b);
1388 namePlot="histDispl_Kcharg_HF_From_b_Bin"; namePlot+=i;
1389 TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1390 hDisplHadr_HF_b->SetMinimum(0);
1391 fOutputStudy->Add(hDisplHadr_HF_b);
1393 namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1394 TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1395 hDisplCharg_b->SetMinimum(0);
1396 fOutputStudy->Add(hDisplCharg_b);
1398 namePlot="histDispl_Charg_HF_From_b_Bin"; namePlot+=i;
1399 TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1400 hDisplCharg_HF_b->SetMinimum(0);
1401 fOutputStudy->Add(hDisplCharg_HF_b);
1403 //origin of tracks histos
1404 namePlot="histOrig_Charg_Bin"; namePlot+=i;
1405 TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1406 hOrigin_Charm->SetMinimum(0);
1407 hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1408 hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1409 hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1410 hOrigin_Charm->GetXaxis()->SetBinLabel(4,"c hadr.");
1411 hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->#");
1412 hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1413 hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->#");
1414 hOrigin_Charm->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1415 hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1416 fOutputStudy->Add(hOrigin_Charm);
1418 namePlot="histOrig_Kcharg_Bin"; namePlot+=i;
1419 TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1420 hOrigin_Kcharg->SetMinimum(0);
1421 hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1422 hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1423 hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1424 hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"c hadr.");
1425 hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->#");
1426 hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1427 hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->#");
1428 hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1429 hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1430 fOutputStudy->Add(hOrigin_Kcharg);
1432 namePlot="histOrig_K0_Bin"; namePlot+=i;
1433 TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1434 hOrigin_K->SetMinimum(0);
1435 hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1436 hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1437 hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1438 hOrigin_K->GetXaxis()->SetBinLabel(4,"c hadr.");
1439 hOrigin_K->GetXaxis()->SetBinLabel(5,"B->#");
1440 hOrigin_K->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1441 hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->#");
1442 hOrigin_K->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1443 hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1444 fOutputStudy->Add(hOrigin_K);
1448 //origin of D0 histos
1449 namePlot="histOrig_D0_Bin"; namePlot+=i;
1450 TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1451 hOrigin_D0->SetMinimum(0);
1452 hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1453 hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1454 fOutputStudy->Add(hOrigin_D0);
1459 //________________________________________________________________________
1460 void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1462 // Method for correlations D0-hadrons study
1464 Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1465 Double_t mD0, mD0bar;
1466 Int_t origD0 = 0, PDGD0 = 0, ptbin = 0;
1467 d->InvMassD0(mD0,mD0bar);
1468 Double_t mInv[2] = {mD0, mD0bar};
1469 ptbin = PtBinCorr(d->Pt());
1471 if(ptbin < 0) return;
1473 //Fill of D0 phi distribution
1474 if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1479 origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1480 PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
1481 switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
1484 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1488 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1495 Double_t highPt = 0; Double_t lead[4] = {0,0,0,1}; //infos for leading particle (pt,deltaphi)
1497 //loop over the tracks in the pool
1498 Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1499 Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1500 Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1502 Int_t NofEventsinPool = 1;
1504 NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1506 AliInfo("Mixed event analysis: track pool is not ready");
1507 NofEventsinPool = 0;
1512 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1513 Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1514 if(!analyzetracksTr) {
1515 AliInfo("AliHFCorrelator::Cannot process the track array");
1519 for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1521 Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1522 if(!runcorrelation) continue;
1524 AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1527 if(!track->CheckSoftPi()) { //removal of soft pions
1528 if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1529 if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1531 } else { //not a soft pion
1532 if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1533 if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1535 Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1536 if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1538 if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1540 Double_t effTr = track->GetWeight(); //extract track efficiency
1541 Double_t effD0 = fD0Eff.at(ptbin);
1542 Double_t eff = effTr*effD0;
1544 FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks, weight = 1./eff
1546 if(!fMixing) N_Charg++;
1548 //retrieving leading info...
1549 if(track->Pt() > highPt) {
1550 if(fReadMC && track->GetLabel()<1) continue;
1551 if(fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) continue;
1552 lead[0] = fCorrelatorTr->GetDeltaPhi();
1553 lead[1] = fCorrelatorTr->GetDeltaEta();
1554 lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1555 lead[3] = 1./track->GetWeight(); //weight is 1./efficiency
1556 highPt = track->Pt();
1559 } // end of tracks loop
1560 } //end of event loop for fCorrelatorTr
1563 NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1565 AliInfo("Mixed event analysis: K+/- pool is not ready");
1566 NofEventsinPool = 0;
1570 //Charged Kaons loop
1571 for (Int_t jMix = 0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1572 Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1573 if(!analyzetracksKc) {
1574 AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1578 for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1580 Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1581 if(!runcorrelation) continue;
1583 AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1586 if(!kCharg->CheckSoftPi()) { //removal of soft pions
1587 if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1588 if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1591 if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1592 if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1594 Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1595 if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1597 if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1599 FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1601 if(!fMixing) N_KCharg++;
1603 } // end of charged kaons loop
1604 } //end of event loop for fCorrelatorKc
1607 NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1609 AliInfo("Mixed event analysis: K0 pool is not ready");
1610 NofEventsinPool = 0;
1615 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1616 Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1617 if(!analyzetracksK0) {
1618 AliInfo("AliHFCorrelator::Cannot process the K0 array");
1622 for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1624 Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1625 if(!runcorrelation) continue;
1627 AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1629 if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1631 FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1633 if(!fMixing) N_Kaons++;
1635 } // end of charged kaons loop
1636 } //end of event loop for fCorrelatorK0
1638 Double_t fillSpLeadD0[3] = {lead[0],mD0,lead[1]};
1639 Double_t fillSpLeadD0bar[3] = {lead[0],mD0bar,lead[1]};
1641 //leading track correlations fill
1644 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1645 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
1646 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
1647 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1648 if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1649 if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //non HF
1651 if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1652 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1653 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
1654 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1655 if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1656 if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //non HF
1659 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);
1660 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1663 //Fill of count histograms
1664 if (!fAlreadyFilled) {
1665 ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1666 ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1667 ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1671 fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1675 //________________________________________________________________________
1676 void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCParticle* d, TClonesArray* mcArray) {
1678 // Method for correlations D0-hadrons study
1680 Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1681 Double_t mD0 = 1.864, mD0bar = 1.864;
1682 Double_t mInv[2] = {mD0, mD0bar};
1683 Int_t origD0 = 0, PDGD0 = 0;
1684 Int_t ptbin = PtBinCorr(d->Pt());
1686 if(ptbin < 0) return;
1688 //Fill of D0 phi distribution
1689 if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1693 origD0=CheckD0Origin(mcArray,d);
1694 PDGD0 = d->GetPdgCode();
1695 switch (CheckD0Origin(mcArray,d)) {
1698 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1702 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1708 Double_t highPt = 0; Double_t lead[3] = {0,0,0}; //infos for leading particle (pt,deltaphi)
1710 //loop over the tracks in the pool
1711 Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1712 Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1713 Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1715 Int_t NofEventsinPool = 1;
1717 NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1719 AliInfo("Mixed event analysis: track pool is not ready");
1720 NofEventsinPool = 0;
1725 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1727 Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1728 if(!analyzetracksTr) {
1729 AliInfo("AliHFCorrelator::Cannot process the track array");
1733 for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1735 Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1736 if(!runcorrelation) continue;
1738 AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1739 if(track->GetLabel()<0) continue;
1740 if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1741 if(track->Pt() < 0.3 || TMath::Abs(track->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1742 if(!fMixing) N_Charg++;
1744 AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1745 if(!trkMC) continue;
1747 if (!trkMC->IsPhysicalPrimary()) continue; //reject material budget, or other fake tracks
1748 if (IsDDaughter(d,trkMC,mcArray)) continue;
1749 FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
1751 //retrieving leading info...
1752 if(track->Pt() > highPt) {
1753 lead[0] = fCorrelatorTr->GetDeltaPhi();
1754 lead[1] = fCorrelatorTr->GetDeltaEta();
1755 lead[2] = fReadMC ? CheckTrackOrigin(mcArray,trkMC) : 0;
1756 highPt = track->Pt();
1759 } // end of tracks loop
1760 } //end of event loop for fCorrelatorTr
1763 NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1765 AliInfo("Mixed event analysis: K+/- pool is not ready");
1766 NofEventsinPool = 0;
1770 //Charged Kaons loop
1771 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1772 Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1773 if(!analyzetracksKc) {
1774 AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1778 for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1780 Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1781 if(!runcorrelation) continue;
1783 AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1784 if(kCharg->GetLabel()<1) continue;
1785 if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1786 if(TMath::Abs(kCharg->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1787 if(!fMixing) N_KCharg++;
1789 AliAODMCParticle *kChargMC = (AliAODMCParticle*)mcArray->At(kCharg->GetLabel());
1790 if(!kChargMC) continue;
1792 if (IsDDaughter(d,kChargMC,mcArray)) continue;
1793 FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1795 } // end of charged kaons loop
1796 } //end of event loop for fCorrelatorKc
1799 NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1801 AliInfo("Mixed event analysis: K0 pool is not ready");
1802 NofEventsinPool = 0;
1807 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1808 Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1809 if(!analyzetracksK0) {
1810 AliInfo("AliHFCorrelator::Cannot process the K0 array");
1814 for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1816 Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1817 if(!runcorrelation) continue;
1819 AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1820 if(k0->GetLabel()<1) continue;
1821 if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1822 if(TMath::Abs(k0->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1824 AliAODMCParticle *k0MC = (AliAODMCParticle*)mcArray->At(k0->GetLabel());
1827 if (IsDDaughter(d,k0MC,mcArray)) continue;
1828 FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1830 if(!fMixing) N_Kaons++;
1832 } // end of charged kaons loop
1833 } //end of event loop for fCorrelatorK0
1835 Double_t fillSpLeadMC[3] = {lead[0],mD0,lead[1]}; //mD0 = mD0bar = 1.864
1837 //leading track correlations fill
1839 if(d->GetPdgCode()==421) { //D0
1840 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
1841 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
1842 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1843 if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1844 if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC); //non HF
1846 if(d->GetPdgCode()==-421) { //D0bar
1847 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
1848 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
1849 if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1850 if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1851 if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC); //non HF
1854 //Fill of count histograms
1855 if (!fAlreadyFilled) {
1856 ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1857 ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1858 ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1861 fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1866 //________________________________________________________________________
1867 void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Double_t mInv[], Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t ptbin, Int_t type, Double_t wg) {
1869 //fills the THnSparse for correlations, calculating the variables
1872 //Initialization of variables
1873 Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
1877 if (fReadMC && track->GetLabel()<1) return;
1878 if (fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) return;
1879 Double_t ptTrack = track->Pt();
1880 Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
1881 Double_t phiTr = track->Phi();
1882 Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1884 TString part = "", orig = "";
1889 deltaphi = fCorrelatorTr->GetDeltaPhi();
1890 deltaeta = fCorrelatorTr->GetDeltaEta();
1895 deltaphi = fCorrelatorKc->GetDeltaPhi();
1896 deltaeta = fCorrelatorKc->GetDeltaEta();
1901 deltaphi = fCorrelatorK0->GetDeltaPhi();
1902 deltaeta = fCorrelatorK0->GetDeltaEta();
1907 if(fMixing == kSE) {
1909 //Fixes limits; needed to include overflow into THnSparse projections!
1910 Double_t pTorig = track->Pt();
1911 Double_t d0orig = track->GetImpPar();
1912 Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1913 Double_t displLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
1914 Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
1915 if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
1916 if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
1917 if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1918 if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1920 //variables for filling histos
1921 Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
1922 Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
1923 Double_t fillSpWeigD0[3] = {deltaphi,mD0,deltaeta};
1924 Double_t fillSpWeigD0bar[3] = {deltaphi,mD0bar,deltaeta};
1927 //sparse fill for data (tracks, K+-, K0) + weighted
1928 if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1929 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1930 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1932 if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1933 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1934 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1936 if(!fAlreadyFilled) {
1937 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1938 ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1944 if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
1946 //sparse fill for data (tracks, K+-, K0) + weighted
1947 if(PdgD0==421) { //D0 (from MCTruth)
1948 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1949 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1950 if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1951 if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1952 if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1953 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1954 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1955 if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1956 if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1957 if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1959 if(PdgD0==-421) { //D0bar
1960 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1961 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1962 if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1963 if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1964 if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1965 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1966 ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1967 if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1968 if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1969 if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1971 if(!fAlreadyFilled) {
1972 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1973 if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
1974 if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
1975 ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1976 ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1977 ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
1978 ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1984 if(fMixing == kME) {
1986 //Fixes limits; needed to include overflow into THnSparse projections!
1987 Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
1988 if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1989 if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1991 //variables for filling histos
1992 Double_t fillSpPhiD0[3] = {deltaphi,mD0,deltaeta};
1993 Double_t fillSpPhiD0bar[3] = {deltaphi,mD0bar,deltaeta};
1996 //sparse fill for data (tracks, K+-, K0)
1997 if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1998 if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2001 //sparse fill for data (tracks, K+-, K0)
2002 if(PdgD0==421) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2003 if(PdgD0==-421) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2011 //_________________________________________________________________________________________________
2012 Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2014 // checks on particle (#) origin:
2018 // 3) c hadronization
2020 // 5) B->X-># (X!=D)
2023 // 8) b hadronization
2025 if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
2027 Int_t pdgGranma = 0;
2029 mother = mcPartCandidate->GetMother();
2031 Int_t abspdgGranma =0;
2032 Bool_t isFromB=kFALSE;
2033 Bool_t isDdaugh=kFALSE;
2034 Bool_t isDchaindaugh=kFALSE;
2035 Bool_t isBdaugh=kFALSE;
2036 Bool_t isBchaindaugh=kFALSE;
2037 Bool_t isQuarkFound=kFALSE;
2039 if (mother<0) return -1;
2040 while (mother >= 0){
2042 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2044 pdgGranma = mcMoth->GetPdgCode();
2045 abspdgGranma = TMath::Abs(pdgGranma);
2046 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2047 isBchaindaugh=kTRUE;
2048 if(istep==1) isBdaugh=kTRUE;
2050 if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
2051 isDchaindaugh=kTRUE;
2052 if(istep==1) isDdaugh=kTRUE;
2054 if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
2055 mother = mcMoth->GetMother();
2057 AliError("Failed casting the mother particle!");
2062 //decides what to return based on the flag status
2064 if(!isFromB) { //charm
2065 if(isDdaugh) return 1; //charm immediate
2066 else if(isDchaindaugh) return 2; //charm chain
2067 else return 3; //charm hadronization
2070 if(isBdaugh) return 4; //b immediate
2071 else if(isBchaindaugh) { //b chain
2073 if(isDdaugh) return 6; //d immediate
2076 else return 5; //b, not d
2078 else return 8; //b hadronization
2081 else if(!isDdaugh && !isDchaindaugh && !isBdaugh && !isBchaindaugh) return 0; //no HF decay
2082 //in this case, it's !isQuarkFound, but not in 100% cases it's a non HF particle!
2083 //rarely you can find a D/B meson which comes from a -1! It isn't a Non-HF, in that case! And I'll return -1...
2085 return -1; //some problem spotted
2087 //________________________________________________________________________
2088 Bool_t AliAnalysisTaskSED0Correlations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
2090 //Daughter removal in MCKine case
2091 Bool_t isDaughter = kFALSE;
2092 Int_t labelD0 = d->GetLabel();
2094 Int_t mother = track->GetMother();
2096 //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
2098 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother)); //it's the mother of the track!
2100 if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
2101 mother = mcMoth->GetMother(); //goes back by one
2103 AliError("Failed casting the mother particle!");
2111 //________________________________________________________________________
2112 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
2114 //give the pt bin where the pt lies.
2117 if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
2120 while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
2125 //---------------------------------------------------------------------------
2126 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
2129 // Selection for K0 hypotheses
2130 // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
2131 // 2 = no previous selections
2133 if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
2135 AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
2136 AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
2138 if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
2139 if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
2142 //This part removes double counting for swapped tracks!
2143 Int_t i = 0; //while loop (until the last-written entry pair of ID!
2144 while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
2145 if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
2146 (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
2149 idArrayV0[i][0]=v0Daug1->GetID();
2150 idArrayV0[i][1]=v0Daug2->GetID();
2155 //________________________________________________________________________
2156 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
2158 cout << "--------------------------\n";
2159 cout << "PtBins = " << fNPtBinsCorr << "\n";
2160 cout << "PtBin limits--------------\n";
2161 for (int i=0; i<fNPtBinsCorr; i++) {
2162 cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
2164 cout << "\n--------------------------\n";
2165 cout << "PtBin tresh. tracks low---\n";
2166 for (int i=0; i<fNPtBinsCorr; i++) {
2167 cout << fPtThreshLow.at(i) << "; ";
2169 cout << "PtBin tresh. tracks up----\n";
2170 for (int i=0; i<fNPtBinsCorr; i++) {
2171 cout << fPtThreshUp.at(i) << "; ";
2173 cout << "D0 efficiencies:----\n";
2174 for (int i=0; i<fNPtBinsCorr; i++) {
2175 cout << fD0Eff.at(i) << "; ";
2177 cout << "\n--------------------------\n";
2178 cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
2179 cout << "--------------------------\n";
2180 cout << "MC Truth = "<<fReadMC<<" - MC Reco: Trk = "<<fRecoTr<<", D0 = "<<fRecoD0<<"\n";
2181 cout << "--------------------------\n";
2182 cout << "Sel of Event tpye (PP,GS,FE,...)= "<<fSelEvType<<"\n";
2183 cout << "--------------------------\n";
2184 cout << "Ev Mixing = "<<fMixing<<"\n";
2185 cout << "--------------------------\n";