]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add p conserv at kine level, flag to speed up the analysis (Fabio C.)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskDStarCorrelations.cxx
CommitLineData
5d3cf93b 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//
16//
17// Base class for DStar - Hadron Correlations Analysis
18//
19//-----------------------------------------------------------------------
20//
21//
22// Author S.Bjelogrlic
23// Utrecht University
24// sandro.bjelogrlic@cern.ch
25//
26//-----------------------------------------------------------------------
27
5d3b1a33 28/* $Id: AliAnalysisTaskDStarCorrelations.cxx 65139 2013-11-25 14:47:45Z fprino $ */
5d3cf93b 29
30//#include <TDatabasePDG.h>
31#include <TParticle.h>
32#include <TVector3.h>
33#include <TChain.h>
34#include "TROOT.h"
35
36#include "AliAnalysisTaskDStarCorrelations.h"
37#include "AliRDHFCutsDStartoKpipi.h"
38#include "AliHFAssociatedTrackCuts.h"
39#include "AliAODRecoDecay.h"
40#include "AliAODRecoCascadeHF.h"
41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliAODPidHF.h"
43#include "AliVParticle.h"
44#include "AliAnalysisManager.h"
45#include "AliAODInputHandler.h"
46#include "AliAODHandler.h"
47#include "AliESDtrack.h"
48#include "AliAODMCParticle.h"
49#include "AliNormalizationCounter.h"
50#include "AliReducedParticle.h"
51#include "AliHFCorrelator.h"
52#include "AliAODMCHeader.h"
53#include "AliEventPoolManager.h"
54#include "AliVertexingHFUtils.h"
55
56using std::cout;
57using std::endl;
58
59
60ClassImp(AliAnalysisTaskDStarCorrelations)
61
62
63//__________________________________________________________________________
64AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
65AliAnalysisTaskSE(),
66fhandler(0x0),
67fmcArray(0x0),
68fCounter(0x0),
69fCorrelator(0x0),
70fselect(0),
71fmontecarlo(kFALSE),
72fmixing(kFALSE),
73fFullmode(kFALSE),
74fSystem(pp),
75fEfficiencyVariable(kNone),
5d3b1a33 76fBkgMethod(kDZeroSB),
5d3cf93b 77fReco(kTRUE),
78fUseEfficiencyCorrection(kFALSE),
79fUseDmesonEfficiencyCorrection(kFALSE),
80fUseCentrality(kFALSE),
81fUseHadronicChannelAtKineLevel(kFALSE),
82fPhiBins(32),
83fEvents(0),
84fDebugLevel(0),
85fDisplacement(0),
86fDim(0),
87fNofPtBins(0),
88fMaxEtaDStar(0.9),
89fDMesonSigmas(0),
90
91fD0Window(0),
92
93fOutput(0x0),
5d3b1a33 94fDmesonOutput(0x0),
95fTracksOutput(0x0),
96fEMOutput(0x0),
97fCorrelationOutput(0x0),
5d3cf93b 98fOutputMC(0x0),
5d3b1a33 99
5d3cf93b 100fCuts(0),
5d3b1a33 101fAssocCuts(0),
5d3cf93b 102fUtils(0),
103fTracklets(0),
104fDeffMapvsPt(0),
105fDeffMapvsPtvsMult(0),
106fDeffMapvsPtvsEta(0)
107{
108 SetDim();
109 // default constructor
110}
111
112//__________________________________________________________________________
113AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :
114AliAnalysisTaskSE(name),
115
116fhandler(0x0),
117fmcArray(0x0),
118fCounter(0x0),
119fCorrelator(0x0),
120fselect(0),
121fmontecarlo(kFALSE),
122fmixing(kFALSE),
123fFullmode(mode),
124fSystem(syst),
125fEfficiencyVariable(kNone),
5d3b1a33 126fBkgMethod(kDZeroSB),
5d3cf93b 127fReco(kTRUE),
128fUseEfficiencyCorrection(kFALSE),
129fUseDmesonEfficiencyCorrection(kFALSE),
130fUseCentrality(kFALSE),
131fUseHadronicChannelAtKineLevel(kFALSE),
132fPhiBins(32),
133fEvents(0),
134fDebugLevel(0),
135fDisplacement(0),
136fDim(0),
137fNofPtBins(0),
138fMaxEtaDStar(0.9),
139fDMesonSigmas(0),
140fD0Window(0),
141
142fOutput(0x0),
5d3b1a33 143fDmesonOutput(0x0),
144fTracksOutput(0x0),
145fEMOutput(0x0),
146fCorrelationOutput(0x0),
5d3cf93b 147fOutputMC(0x0),
5d3b1a33 148
5d3cf93b 149fCuts(0),
5d3b1a33 150fAssocCuts(AsscCuts),
5d3cf93b 151fUtils(0),
152fTracklets(0),
153fDeffMapvsPt(0),
154fDeffMapvsPtvsMult(0),
155fDeffMapvsPtvsEta(0)
156{
157 Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
158
159 SetDim();
5d3b1a33 160 if(fSystem == pp) fUseCentrality = kFALSE; else fUseCentrality = kTRUE;
161
162 if(fSystem == AA) fBkgMethod = kDStarSB; else fBkgMethod = kDZeroSB;
5d3cf93b 163
164 fCuts=cuts;
165 fNofPtBins= fCuts->GetNPtBins();
166 //cout << "Enlarging the DZero window " << endl;
167 EnlargeDZeroMassWindow();
168 // cout << "Done" << endl;
169
170
171 DefineInput(0, TChain::Class());
172
173 DefineOutput(1,TList::Class()); // histos from data and MC
174 DefineOutput(2,TList::Class()); // histos from MC
5d3b1a33 175 DefineOutput(3,TList::Class()); // histos from data and MC
176 DefineOutput(4,TList::Class()); // histos from MC
177 DefineOutput(5,TList::Class()); // histos from data and MC
178 DefineOutput(6,TList::Class()); // histos from data and MC
179 DefineOutput(7,AliNormalizationCounter::Class()); // normalization
180
181 DefineOutput(8,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
182 DefineOutput(9,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
5d3cf93b 183}
184
185//__________________________________________________________________________
186
187AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
188 //
189 // destructor
190 //
191
192 Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");
193
194 if(fhandler) {delete fhandler; fhandler = 0;}
195 //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}
196 if(fmcArray) {delete fmcArray; fmcArray = 0;}
197 if(fCounter) {delete fCounter; fCounter = 0;}
198 if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
199 if(fOutput) {delete fOutput; fOutput = 0;}
200 if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
201 if(fCuts) {delete fCuts; fCuts = 0;}
5d3b1a33 202 if(fAssocCuts) {delete fAssocCuts; fAssocCuts=0;}
5d3cf93b 203 if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}
204 if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}
205 if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}
206
207}
208
209//___________________________________________________________
210void AliAnalysisTaskDStarCorrelations::Init(){
211 //
212 // Initialization
213 //
214 if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
215
216 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
217
218
219
220 // Post the D* cuts
5d3b1a33 221 PostData(8,copyfCuts);
5d3cf93b 222
223 // Post the hadron cuts
5d3b1a33 224 PostData(9,fAssocCuts);
5d3cf93b 225
226 return;
227}
228
229
230//_________________________________________________
231void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
232 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
233
234 //slot #1
235 //OpenFile(0);
236 fOutput = new TList();
237 fOutput->SetOwner();
238
239 fOutputMC = new TList();
240 fOutputMC->SetOwner();
5d3b1a33 241
242 fDmesonOutput = new TList();
243 fDmesonOutput->SetOwner();
244
245 fTracksOutput = new TList();
246 fTracksOutput->SetOwner();
247
248 fEMOutput = new TList();
249 fEMOutput->SetOwner();
250
251 fCorrelationOutput = new TList();
252 fCorrelationOutput->SetOwner();
5d3cf93b 253
254 // define histograms
255 DefineHistoForAnalysis();
256 DefineThNSparseForAnalysis();
257
258
5d3b1a33 259
260
261 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(6)->GetContainer()->GetName()));
5d3cf93b 262 fCounter->Init();
263
264 Double_t Pi = TMath::Pi();
5d3b1a33 265 //AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t useCentrality, AliRDHFCuts *cutObject)
266 fCorrelator = new AliHFCorrelator("Correlator",fAssocCuts,fUseCentrality,fCuts); // fAssocCuts is the hadron cut object, fSystem to switch between pp or PbPb
5d3cf93b 267 fCorrelator->SetDeltaPhiInterval( -0.5*Pi, 1.5*Pi); // set correct phi interval
268 //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval
269 fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
270 fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
271 fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
272 fCorrelator->SetUseMC(fmontecarlo);
273 fCorrelator->SetUseReco(fReco);
274 // fCorrelator->SetKinkRemoval(kTRUE);
275 Bool_t pooldef = fCorrelator->DefineEventPool();
276
277 if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
278
279 fUtils = new AliAnalysisUtils();
280
281
282
283 PostData(1,fOutput); // set the outputs
5d3b1a33 284 PostData(2,fDmesonOutput); // set the outputs
285 PostData(3,fTracksOutput); // set the outputs
286 PostData(4,fEMOutput); // set the outputs
287 PostData(5,fCorrelationOutput); // set the outputs
288 PostData(6,fOutputMC); // set the outputs
289 PostData(7,fCounter); // set the outputs
5d3cf93b 290}
291
5d3b1a33 292//________________________________________ user exec ____________
5d3cf93b 293void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
5d3b1a33 294 // cout << "Task debug check 1 " << endl;
295 // ********************************************** LOAD THE EVENT ****************************************************
5d3cf93b 296
5d3b1a33 297 if (!fInputEvent) {
298 Error("UserExec","NO EVENT FOUND!");
5d3cf93b 299 return;
5d3b1a33 300 }
301
302 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
303 if(!aodEvent){
304 AliError("AOD event not found!");
305 return;
306 }
5d3cf93b 307
308 fEvents++; // event counter
309 ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
310
5d3b1a33 311 fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); // store event in AliNormalizationCounter
5d3cf93b 312
313 // load MC array
314 fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
315 if(fmontecarlo && !fmcArray){
5d3b1a33 316 AliError("Array of MC particles not found");
317 return;
5d3cf93b 318 }
319
5d3cf93b 320 // ********************************************** START EVENT SELECTION ****************************************************
321
322 Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
323
324 if(!isEvSel) {
5d3b1a33 325
326 if(fCuts->IsEventRejectedDueToPileup()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2); cout << "Reject PILEUP" << endl;}
327 if(fCuts->IsEventRejectedDueToCentrality()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3); cout << "Reject CENTRALITY" << endl;}
328 if(fCuts->IsEventRejectedDueToNotRecoVertex()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4); cout << "Reject NOT RECO VTX" << endl;}
329 if(fCuts->IsEventRejectedDueToVertexContributors()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5); cout << "Reject VTX CONTRIB" << endl;}
330 if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6); cout << "Reject VTX no fid reg " << endl;}
331 if(fCuts->IsEventRejectedDueToTrigger()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7); cout << "Reject TRIGGER" << endl;}
332 if(fCuts->IsEventRejectedDuePhysicsSelection()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8); cout << "Reject PHYS SEL" << endl;}
333
334 return;
5d3cf93b 335 }
5d3cf93b 336 // added event selection for pA
337
338 if(fSystem == pA){
5d3b1a33 339
340 if(fUtils->IsFirstEventInChunk(aodEvent)) {
341 AliInfo("Rejecting the event - first in the chunk");
342 ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);
343 return;
344 }
345 if(!fUtils->IsVertexSelected2013pA(aodEvent)) {
346 ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);
347 AliInfo("Rejecting the event - bad vertex");
348 return;
5d3cf93b 349 }
5d3cf93b 350 }
5d3b1a33 351
352 fCorrelator->SetAODEvent(aodEvent); // set the event in the correlator class
5d3cf93b 353
5d3b1a33 354 Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event within the pool definition
355 if(!correlatorON) {
356 AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
357 return; // if not the case, skips the event
358 }
5d3cf93b 359
5d3b1a33 360 ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); // count events that are passing the event selection
5d3cf93b 361
5d3cf93b 362
5d3cf93b 363
5d3b1a33 364 // cout << "Task debug check 2 " << endl;
365 // ********************************************** CENTRALITY/MULTIPLICITY ****************************************************
366
367
368 Double_t MultipOrCent = fCorrelator->GetCentrality(); // returns centrality or multiplicity
5d3cf93b 369
5d3b1a33 370
371
372 // ********************************************** MC SETUP ****************************************************
373
374 if(fmontecarlo) {
375
376 fCorrelator->SetMCArray(fmcArray); // export MC array into correlatior class
377
378 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
379 if (fmontecarlo && !mcHeader) {
380 AliError("Could not find MC Header in AOD");
381 return;
382 }
383
384 // select MC event type (PP, GS etc) - those are defined in the associated tracks cut object
385 Bool_t isMCeventgood = kFALSE;
386
387
388 Int_t eventType = mcHeader->GetEventType();
389 Int_t NMCevents = fAssocCuts->GetNofMCEventType();
390
391 for(Int_t k=0; k<NMCevents; k++){
392 Int_t * MCEventType = fAssocCuts->GetMCEventType();
393
394 if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
395 ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
396 }
397
398 if(NMCevents && !isMCeventgood){
5d3cf93b 399 if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
5d3b1a33 400 return;
401 }
402
403
404 }// end if fmontecarlo
405
406
407 // ********************************************** EVENT PROPERTIES ****************************************************
408 // checks on vertex and multiplicity of the event
5d3cf93b 409 AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
410 Double_t zVtxPosition = vtx->GetZ(); // zvertex
411
5d3b1a33 412
413
414 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
415
416 if(!fmixing) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);
417 if(!fmixing) ((TH1D*)fOutput->FindObject("MultiplicityOnlyCheck"))->Fill(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
418
419
420 // Int_t poolbin = fAssocCuts->GetPoolBin(MultipOrCent, zVtxPosition); // get the event pool bin - commented for the moment - to be checked
421
422
423
424 // ********************************************** D * selection ****************************************************
425
426
427 TClonesArray *arrayDStartoD0pi=0;
5d3cf93b 428 if(!aodEvent && AODEvent() && IsStandardAOD()) {
5d3b1a33 429 // In case there is an AOD handler writing a standard AOD, use the AOD
430 // event in memory rather than the input (ESD) event.
431 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
432 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
433 // have to taken from the AOD event hold by the AliAODExtension
434 AliAODHandler* aodHandler = (AliAODHandler*)
5d3cf93b 435 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
5d3b1a33 436 if(aodHandler->GetExtensions()) {
437 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
438 AliAODEvent *aodFromExt = ext->GetAOD();
439 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
440 }
5d3cf93b 441 } else {
5d3b1a33 442 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
5d3cf93b 443 }
5d3cf93b 444
5d3cf93b 445
5d3cf93b 446
5d3b1a33 447 // D* related variables
5d3cf93b 448
5d3b1a33 449 Double_t ptDStar; // D* pt
450 Double_t phiDStar; // D* phi
451 Double_t etaDStar; // D* eta
452 Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; // flags for selection
453 Double_t invMassDZero; // D0 inv mass
454 Double_t deltainvMDStar; // D* - D0 invarian mass
455
456 Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
5d3cf93b 457 Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
458
459
460 //MC tagging for DStar
461 //D* and D0 prongs needed to MatchToMC method
462 Int_t pdgDgDStartoD0pi[2]={421,211};
463 Int_t pdgDgD0toKpi[2]={321,211};
464
5d3b1a33 465 //Bool_t isDStarCand = kFALSE;
5d3cf93b 466 Bool_t isDfromB = kFALSE;
467 Bool_t isEventMixingFilledPeak = kFALSE;
468 Bool_t isEventMixingFilledSB = kFALSE;
469 Bool_t EventHasDStarCandidate = kFALSE;
470 Bool_t EventHasDZeroSideBandCandidate = kFALSE;
471 Bool_t EventHasDStarSideBandCandidate = kFALSE;
5d3b1a33 472 Bool_t isTrackForPeakFilled = kFALSE;
473 Bool_t isTrackForSBFilled = kFALSE;
5d3cf93b 474 //loop on D* candidates
475
476 Int_t looponDCands = 0;
5d3b1a33 477 if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); // load N of D* candidates from AOD array
478 if(!fReco) looponDCands = fmcArray->GetEntriesFast(); // load N of D* candidates from MC array
5d3cf93b 479
5d3b1a33 480 Int_t nOfDStarCandidates = 0; // D candidates counter
481 Int_t nOfSBCandidates = 0; // SB candidates counter
5d3cf93b 482
5d3b1a33 483 Double_t DmesonEfficiency = 1.; // Efficiency of D meson
484 Double_t DmesonWeight = 1.; // Applyed weight
485 Double_t efficiencyvariable = -999; // Variable to be used (defined by the AddTask)
486 Double_t dmDStarWindow = 0.0019/3; // DStar window
5d3cf93b 487
5d3b1a33 488 // cout << "Task debug check 3 " << endl;
489 // loop over D meson candidates
490 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
491
492 // initialize variables
493 isInPeak = kFALSE;
494 isInDStarSideBand = kFALSE;
495 isInDZeroSideBand = kFALSE;
496
497 EventHasDStarCandidate = kFALSE;
498 EventHasDZeroSideBandCandidate = kFALSE;
499 EventHasDStarSideBandCandidate = kFALSE;
500
501
502 isDStarMCtag = kFALSE;
503 isDfromB = kFALSE;
504 ptDStar = -123.4;
505 phiDStar = -999;
506 etaDStar = -56.;
507 invMassDZero = - 999;
508 deltainvMDStar = -998;
509 AliAODRecoCascadeHF* dstarD0pi;
510 AliAODRecoDecayHF2Prong* theD0particle;
511 AliAODMCParticle* DStarMC=0x0;
512 Short_t daughtercharge = -2;
513 Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
514 Int_t trackiddaugh1 = -1;
515 Int_t trackidsoftPi = -1;
516 Int_t ptbin = -1;
517
518 // ============================================ using reconstruction on Data or MC
519 if(fReco){
520 // cout << "Task debug check 4 " << endl;
521 // get the D objects
522 dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
523 if(!dstarD0pi->GetSecondaryVtx()) continue;
524 theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
525 if (!theD0particle) continue;
5d3cf93b 526
5d3cf93b 527
5d3b1a33 528 // apply topological cuts
5d3cf93b 529
5d3b1a33 530 // track quality cuts
531 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
532 // region of interest + topological cuts + PID
533 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
5d3cf93b 534
5d3b1a33 535 //apply topological cuts
536 if(!isTkSelected) continue;
537 if(!isSelected) continue;
538 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
539
540 // get D candidate variables
541 ptDStar = dstarD0pi->Pt();
542 phiDStar = dstarD0pi->Phi();
543 etaDStar = dstarD0pi->Eta();
544 if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
545 if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr) efficiencyvariable = MultipOrCent;
546 if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
547 if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
548 if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
549
550
551 if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
5d3cf93b 552
5d3b1a33 553 phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
554 ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
5d3cf93b 555
5d3b1a33 556 cout << "DStar pt = " << ptDStar << endl;
557 cout << "pt bin = " << ptbin << endl;
558 if(ptbin<1) continue;
5d3cf93b 559
5d3b1a33 560 Double_t mD0Window= fD0Window[ptbin]/3;
561 invMassDZero = dstarD0pi->InvMassD0();
562 deltainvMDStar = dstarD0pi->DeltaInvMass();
563
564
565 // get the D meson efficiency
566 DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
567
568 // check this!
569 if(fUseDmesonEfficiencyCorrection){
570 if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
571 else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI
572 if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value
573 else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low
574 }
575 }
576 else DmesonWeight = 1.;
5d3cf93b 577
5d3b1a33 578
579 // using montecarlo on reconstruction
580 Int_t mcLabelDStar = -999;
581 if(fmontecarlo){
582 // find associated MC particle for D* ->D0toKpi
583 mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);
584 if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
585 if(!isDStarMCtag) continue;
586 AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
587 //check if DStar from B
588 Int_t labelMother = MCDStar->GetMother();
589 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
590 if(!mother) continue;
591 Int_t motherPDG =TMath::Abs(mother->PdgCode());
592 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
593 }
5d3cf93b 594
5d3cf93b 595
5d3b1a33 596 // fill mass histograms
5d3cf93b 597
5d3b1a33 598 // cout << "Task debug check 5 " << endl;
599 // fill D0 invariant mass
600 if(!fmixing) {
601 cout << " histo name = " << Form("histDZeroMass_%d",ptbin) << endl;
602 ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
603 // cout << "Task debug check 5.1 " << endl;
604 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
605 } // end if !mixing
5d3cf93b 606
5d3b1a33 607
5d3cf93b 608
5d3cf93b 609
5d3cf93b 610
5d3b1a33 611 // good D0 canidates
612 if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
613 // fill D* invariant mass
614 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
615 // fill D* invariant mass if weighting
616 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
617 isInPeak = kTRUE;
618 // fill other good candidate distributions - pt, phi eta
619 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
620 ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
621 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
622 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
623 nOfDStarCandidates++;
624 EventHasDStarCandidate=kTRUE;
625 } // end if in good D* mass window
626
627 // count D* sideband candidates
628 if(fBkgMethod == kDStarSB ){
629 if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
630 ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
631 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
632 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
633 nOfSBCandidates++;
634 EventHasDStarSideBandCandidate = kTRUE;
635 }
636
637 } // end if using DStar sidebans
638
639
640 }// end good D0
641
642 // cout << "Task debug check 6 " << endl;
643 // Sideband D0
644 if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
645 isInDZeroSideBand = kTRUE;
646 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
647 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
648
649 if(fBkgMethod == kDZeroSB){
650 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
651
652 ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
653 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
654 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
655 nOfSBCandidates++;
656 EventHasDZeroSideBandCandidate = kTRUE;
657 }
658 }
659 }
660 }// end SideBand D0
661 // cout << "Task debug check 7 " << endl;
5d3cf93b 662
5d3b1a33 663 if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
664 if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.151) continue; // skip all D* canidates outside the interesting mass ranges
665 // cout << "Good D* candidate" << endl;
5d3cf93b 666
5d3b1a33 667 // get charge of soft pion
668 daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
5d3cf93b 669
5d3b1a33 670 // get ID of daughters used for reconstruction
671 trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
672 trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
673 trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
674 }// end if reconstruction
675
676
677
678 // ============================================ using MC kinematics only
679 Int_t DStarLabel = -1;
680
681 if(!fReco){ // use pure MC information
5d3cf93b 682
5d3b1a33 683 // get the DStar Particle
684 DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
685 if (!DStarMC) {
686 AliWarning("Careful: DStar MC Particle not found in tree, skipping");
687 continue;
688 }
689 DStarLabel = DStarMC->GetLabel();
690 if(DStarLabel>0)isDStarMCtag = kTRUE;
5d3cf93b 691
5d3b1a33 692 Int_t PDG =TMath::Abs(DStarMC->PdgCode());
693 if(PDG !=413) continue; // skip if it is not a DStar
5d3cf93b 694
5d3b1a33 695 // check fiducial acceptance
696 if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
697 ptbin=fCuts->PtBin(DStarMC->Pt());
5d3cf93b 698
699 //check if DStar from B
700 Int_t labelMother = DStarMC->GetMother();
701 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
5d3b1a33 702 if(!mother) continue;
5d3cf93b 703 Int_t motherPDG =TMath::Abs(mother->PdgCode());
704 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
5d3b1a33 705
706
707 Bool_t isDZero = kFALSE;
708 Bool_t isSoftPi = kFALSE;
709
710 if(fUseHadronicChannelAtKineLevel){
711 //check decay channel on MC ************************************************
5d3cf93b 712 Int_t NDaugh = DStarMC->GetNDaughters();
713 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
714
5d3b1a33 715 for(Int_t i=0; i<NDaugh;i++){ // loop on DStar daughters
716 Int_t daugh_label = DStarMC->GetDaughter(i);
717 AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
718 if(!mcDaughter) continue;
719 Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
720 if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
721
722 if(daugh_pdg == 421) {
723 Int_t NDaughD0 = mcDaughter->GetNDaughters();
724 if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
725 trackiddaugh0 = mcDaughter->GetDaughter(0);
726 trackiddaugh1 = mcDaughter->GetDaughter(1);
727 Bool_t isKaon = kFALSE;
728 Bool_t isPion = kFALSE;
729
730 for(Int_t k=0;k<NDaughD0;k++){
731 Int_t labelD0daugh = mcDaughter->GetDaughter(k);
732 AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
733 if(!mcGrandDaughter) continue;
734 Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
735 if(granddaugh_pdg==321) isKaon = kTRUE;
736 if(granddaugh_pdg==211) isPion = kTRUE;
5d3cf93b 737 }
5d3b1a33 738 if(!isKaon || !isPion) break; // skip if not correct decay channel of D0
739 isDZero = kTRUE;
740 } // end check on Dzero
741
742 if(daugh_pdg == 211) {
743 isSoftPi = kTRUE;
744 daughtercharge = mcDaughter->Charge();
745 trackidsoftPi = daugh_label;}
746 } // end loop on D* daughters
747 if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
748 } // end check decay channel
749
750 ptDStar = DStarMC->Pt();
751 phiDStar = DStarMC->Phi();
752 etaDStar = DStarMC->Eta();
753
754 if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
755
756
757 }// end if pure MC information
758
759
760
761
762
5d3cf93b 763 // getting the number of triggers in the MCtag D* case
5d3b1a33 764 if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar);
765 if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
766 if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
767
768
769 fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
770 fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
771
772
773
774
775 // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
776
777 cout << "Correlating " << endl;
778
779 Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
780
781 // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
782 if(fmixing && !execPool) {
783 AliInfo("Mixed event analysis: pool is not ready");
5d3cf93b 784 if(!isEventMixingFilledPeak && isInPeak) {
5d3b1a33 785 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
786 isEventMixingFilledPeak = kTRUE;
5d3cf93b 787 }
788 if (!isEventMixingFilledSB && isInDZeroSideBand) {
5d3b1a33 789 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
790 isEventMixingFilledSB=kTRUE;
5d3cf93b 791 }
5d3b1a33 792 continue;
793 } // end if pool not ok
794 // if analysis is on mixed event and pool settings are satisfied, fill relevant histograms and continue
795 if(fmixing&&execPool){
5d3cf93b 796 // pool is ready - run checks on bins filling
797 if(!isEventMixingFilledPeak && isInPeak) {
5d3b1a33 798 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
799 if(fFullmode) EventMixingChecks(aodEvent);
800 isEventMixingFilledPeak = kTRUE;
5d3cf93b 801 }
802
803 if(!isEventMixingFilledSB && isInDZeroSideBand) {
5d3b1a33 804 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
805 isEventMixingFilledSB=kTRUE;
5d3cf93b 806 }
5d3b1a33 807 } // end if pool ok
808
809
810
811
812 Int_t NofEventsinPool = 1;
813 if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
814
815 //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
816
817 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
5d3cf93b 818
5d3b1a33 819 Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); // process the associated tracks
820 if(!analyzetracks) {
821 AliInfo("AliHFCorrelator::Cannot process the track array");
822 continue;
823 }
5d3cf93b 824
5d3b1a33 825 Double_t arraytofill[4];
826 // Double_t MCarraytofill[4]; // think about this
5d3cf93b 827 Double_t weight;
828
5d3b1a33 829 Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
5d3cf93b 830
5d3b1a33 831 //************************************************** LOOP ON TRACKS *****************************************************
5d3cf93b 832
5d3b1a33 833 for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
834
835 Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
836 if(!runcorrelation) continue;
837
838 Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
839 Double_t DeltaEta = fCorrelator->GetDeltaEta();
840
841 AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
842 if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
843
844 Double_t ptHad = hadron->Pt();
845 Double_t phiHad = hadron->Phi();
846 phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
847 Double_t etaHad = hadron->Eta();
848 Int_t label = hadron->GetLabel();
849 Int_t trackid = hadron->GetID();
850 Double_t efficiency = hadron->GetWeight();
851
852
853
854
855 if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
856
857 ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
858 ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
859 isTrackForPeakFilled = kTRUE; // makes sure you do not fill twice in case of more candidates
860 }
861
862 if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
863 ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
864 ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
865 isTrackForSBFilled = kTRUE;
866 }
867
868 if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
869 ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
870 ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
871 isTrackForSBFilled = kTRUE;
872 }
873
874
875 weight = 1;
876 if(fUseEfficiencyCorrection && efficiency){
877 weight = DmesonWeight * (1./efficiency);
878 }
879
880
881 arraytofill[0] = DeltaPhi;
882 arraytofill[1] = deltainvMDStar;
883 arraytofill[2] = DeltaEta;
884 arraytofill[3] = ptHad;
885 // arraytofill[4] = poolbin;
886
887
888 // skip the D daughters in the correlation
889 // Bool_t isDdaughter = kFALSE;
890 if(!fmixing && fReco){ // for reconstruction
891 if(trackid == trackiddaugh0) continue;
892 if(trackid == trackiddaugh1) continue;
893 if(trackid == trackidsoftPi) continue;
894 }
895
896 if(!fmixing && !fReco){ // for kinematic MC
897 AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
898 if(!part) continue;
899 if(IsDDaughter(DStarMC, part)) continue;
900 cout << "Skipping DStar daugheter " << endl;
901 }
902
903
904
905 // ============================================= FILL CORRELATION THNSparses ===============================
906
907 // filling signal
908 if(isInPeak){
909 cout << "Filling signal " << endl;
910 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
911 //cout ("CorrelationsDStarHadron_%d",ptbin)
912 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
913 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
914 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
915 }
916
917 // filling bkg
918 if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
919 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
920 cout << "Filling bkg from D* sidebands " << endl;
921 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
922 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
923 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
924
925 } // end if bkg from DStar
926
927 if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
928 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
929 cout << "Filling bkg from Dzero sidebands " << endl;
930 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
931 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
932 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
933
934 } // end if bkg from DZero
935
936
937
938 } // end loop on associated tracks
5d3cf93b 939
5d3b1a33 940 } // end loop on events in event pool
5d3cf93b 941
5d3b1a33 942
943 } // end loop on D mesons
5d3cf93b 944
5d3cf93b 945
5d3cf93b 946
947
5d3b1a33 948 Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
949
950
951 if(!updated) AliInfo("Pool was not updated");
952
5d3cf93b 953
954} //end the exec
955
956//________________________________________ terminate ___________________________
957void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
958{
959 // The Terminate() function is the last function to be called during
960 // a query. It always runs on the client, it can be used to present
961 // the results graphically or save the results to file.
962
963 AliAnalysisTaskSE::Terminate();
964
965 fOutput = dynamic_cast<TList*> (GetOutputData(1));
966 if (!fOutput) {
967 printf("ERROR: fOutput not available\n");
968 return;
969 }
970
971 return;
972}
973//_____________________________________________________
974Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
975
976 //Daughter removal in MCKine case
977 Bool_t isDaughter = kFALSE;
978 Int_t labelD0 = d->GetLabel();
979
980 Int_t mother = track->GetMother();
981
982 //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
983 while (mother > 0){
984 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
985 if (mcMoth){
986 if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
987 mother = mcMoth->GetMother(); //goes back by one
988 } else{
989 AliError("Failed casting the mother particle!");
990 break;
991 }
992 }
993
994 return isDaughter;
995}
996
997//_____________________________________________________
998void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
999
5d3cf93b 1000 Double_t Pi = TMath::Pi();
1001 Int_t nbinscorr = fPhiBins;
5d3b1a33 1002 Double_t lowcorrbin = -0.5*Pi ;
1003 Double_t upcorrbin = 1.5*Pi ;
5d3cf93b 1004
5d3cf93b 1005
5d3b1a33 1006 // create ThNSparses
5d3cf93b 1007
5d3b1a33 1008 Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
5d3cf93b 1009
1010
5d3b1a33 1011 //sparse bins
1012
1013 //1 delta_phi
1014 //2 invariant mass D *
1015 //3 delta eta
1016 //4 track pt
5d3cf93b 1017
5d3cf93b 1018
5d3b1a33 1019 //Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
5d3cf93b 1020
5d3cf93b 1021
5d3b1a33 1022 // for reconstruction on Data and MC
1023 Int_t nbinsSparse[4]= {nbinscorr , 32 , 32, 10};
1024 Double_t binLowLimitSparse[4]={lowcorrbin,0.14314 ,-1.6, 0};
1025 Double_t binUpLimitSparse[4]= {upcorrbin ,0.14794 , 1.6, 5};
5d3cf93b 1026
5d3b1a33 1027 Int_t nbinsSparseDStarSB[4]= {nbinscorr , 40 , 32, 10};
1028 Double_t binLowLimitSparseDStarSB[4]={lowcorrbin,0.14788 ,-1.6, 0};
1029 Double_t binUpLimitSparseDStarSB[4]= {upcorrbin ,0.1504 , 1.6, 5};
5d3cf93b 1030
5d3b1a33 1031 TString signalSparseName = "";
1032 TString bkgSparseName = "";
5d3cf93b 1033
5d3b1a33 1034 THnSparseF * CorrelationsSignal = NULL;
1035 THnSparseF * CorrelationsBkg = NULL;
5d3cf93b 1036
5d3cf93b 1037
5d3b1a33 1038 Float_t * ptbinlims = fCuts->GetPtBinLimits();
1039
1040
5d3cf93b 1041
5d3cf93b 1042
5d3b1a33 1043 for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1044
1045
1046
1047 if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1048
1049
1050
1051 signalSparseName = "CorrelationsDStar";
1052 if(fselect==1) signalSparseName += "Hadron_";
1053 if(fselect==2) signalSparseName += "Kaon_";
1054 if(fselect==3) signalSparseName += "KZero_";
1055
1056 bkgSparseName = "CorrelationsBkg";
1057 if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
1058 if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
1059 if(fselect==1) bkgSparseName += "Hadron_";
1060 if(fselect==2) bkgSparseName += "Kaon_";
1061 if(fselect==3) bkgSparseName += "KZero_";
1062
1063 signalSparseName+=Form("%d",iBin);
1064 bkgSparseName+=Form("%d",iBin);
1065 cout << "ThNSparses name = " << signalSparseName << endl;
1066
1067 // define thnsparses for signal candidates
1068 CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1069 CorrelationsSignal->Sumw2();
1070 fCorrelationOutput->Add(CorrelationsSignal);
1071
1072 // define thnsparses for bkg candidates from DStar
1073 if(fBkgMethod == kDStarSB ){
1074 CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",4,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
1075 CorrelationsBkg->Sumw2();
1076 fCorrelationOutput->Add(CorrelationsBkg);
1077 }
1078
1079 // define thnsparses for bkg candidates from DZero
1080 if(fBkgMethod == kDZeroSB ){
1081 CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1082 CorrelationsBkg->Sumw2();
1083 fCorrelationOutput->Add(CorrelationsBkg);
1084 }
1085
1086 } // end loop on bins
5d3cf93b 1087
5d3cf93b 1088
1089
1090}
5d3b1a33 1091
5d3cf93b 1092//__________________________________________________________________________________________________
1093void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
5d3b1a33 1094
1095 Double_t Pi = TMath::Pi();
5d3cf93b 1096 Int_t nbinscorr = fPhiBins;
5d3b1a33 1097 Double_t lowcorrbin = -0.5*Pi ;
1098 Double_t upcorrbin = 1.5*Pi ;
1099
1100 // event counter
1101 TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
5d3cf93b 1102 NofEvents->GetXaxis()->SetBinLabel(1," All events");
1103 NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1104 NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1105 NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1106 NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1107 NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1108 NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1109 NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1110 NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1111 NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1112 NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1113 fOutput->Add(NofEvents);
5d3cf93b 1114
5d3b1a33 1115 //event properties
1116 TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1117 fOutput->Add(EventPropsCheck);
5d3cf93b 1118
5d3b1a33 1119 //event properties
1120 TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",1000,0,1000);
1121 fOutput->Add(SPDMultiplicty);
5d3cf93b 1122
5d3cf93b 1123
5d3b1a33 1124 // =================================================== D* histograms
1125 TH1F * D0mass = NULL;
1126 TH1F * DStarMass = NULL;
1127 TH1F * DStarFromSBMass = NULL;
5d3cf93b 1128
5d3b1a33 1129 TH1F * D0massWeighted = NULL;
1130 TH1F * DStarMassWeighted = NULL;
1131 TH1F * DStarFromSBMassWeighted = NULL;
5d3cf93b 1132
5d3b1a33 1133 TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "";
5d3cf93b 1134
5d3b1a33 1135 Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1136 Float_t * ptbinlims = fCuts->GetPtBinLimits();
5d3cf93b 1137
5d3b1a33 1138 //GetMinPtCandidate()
1139
5d3cf93b 1140
5d3b1a33 1141 for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1142
1143
1144
1145 if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1146
1147
1148 std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
1149
1150 nameDZeroMass = "histDZeroMass_";
1151 nameDStarMass = "histDStarMass_";
1152 nameDStarFromSBMass = "histDStarFromSBMass_";
1153
1154 nameDZeroMass+=Form("%d",iBin);
1155 nameDStarMass+=Form("%d",iBin);
1156 nameDStarFromSBMass+=Form("%d",iBin);
1157
1158 cout << "D zero histogram: " << nameDZeroMass << endl;
1159
1160 D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1161 DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1162 DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1163
1164 if(!fmixing){
1165 fDmesonOutput->Add(D0mass);
1166 fDmesonOutput->Add(DStarMass);
1167 fDmesonOutput->Add(DStarFromSBMass);
1168 }
1169
1170 // if using D meson efficiency, define weighted histos
1171 if(fUseDmesonEfficiencyCorrection){
1172
1173 nameDZeroMass = "histDZeroMassWeight_";
1174 nameDStarMass = "histDStarMassWeight_";
1175 nameDStarFromSBMass = "histDStarFromSBMassWeight_";
1176
1177 nameDZeroMass+=Form("%d",iBin);
1178 nameDStarMass+=Form("%d",iBin);
1179 nameDStarFromSBMass+=Form("%d",iBin);
1180
1181 D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1182 DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1183 DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1184
1185 if(!fmixing){
1186 fDmesonOutput->Add(D0massWeighted);
1187 fDmesonOutput->Add(DStarMassWeighted);
1188 fDmesonOutput->Add(DStarFromSBMassWeighted);
1189 }
1190 }
1191 }// end loop on pt bins
5d3cf93b 1192
5d3cf93b 1193
5d3b1a33 1194 // pt distributions
1195 TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
1196 fDmesonOutput->Add(RecoPtDStar);
1197 TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
1198 fDmesonOutput->Add(RecoPtBkg);
5d3cf93b 1199
5d3b1a33 1200 // phi distribution
1201 TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1202 TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
5d3cf93b 1203
5d3b1a33 1204 // eta distribution
1205 TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1206 TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
5d3cf93b 1207
5d3b1a33 1208 if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
1209 if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
1210 if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
1211 if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
5d3cf93b 1212
5d3cf93b 1213
5d3b1a33 1214 // single track related histograms
1215 // phi distribution
1216 TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1217 TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
5d3cf93b 1218
5d3b1a33 1219 // eta distribution
1220 TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1221 TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
5d3cf93b 1222
5d3b1a33 1223 if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
1224 if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
1225 if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
1226 if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
1227
5d3cf93b 1228
5d3b1a33 1229 // Montecarlo for D*
1230 TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1231 if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
5d3cf93b 1232
5d3b1a33 1233 TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1234 if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
5d3cf93b 1235
5d3b1a33 1236 TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1237 if(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
5d3cf93b 1238
5d3cf93b 1239
5d3b1a33 1240 // event mixing histograms
5d3cf93b 1241 TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1242 CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1243 CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1244 CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1245 CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1246
5d3b1a33 1247 if(fmixing) fEMOutput->Add(CheckPoolReadiness);
5d3cf93b 1248
5d3b1a33 1249
1250 /* Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
1251 Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
1252 Int_t nPoolBins = NofCentBins*NofZVrtxBins;
5d3cf93b 1253
5d3b1a33 1254
1255 TH1D * PoolBinDistribution = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1256 fEMOutput->Add(PoolBinDistribution);
1257
1258 TH2D * EventDistributionPerPoolBin = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1259 fEMOutput->Add(EventDistributionPerPoolBin);
1260 */
5d3cf93b 1261}
1262
5d3b1a33 1263
5d3cf93b 1264//__________________________________________________________________________________________________
1265void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1266
1267
1268 //Float_t* ptbins = fCuts->GetPtBinLimits();
1269 if(fD0Window) delete fD0Window;
1270 fD0Window = new Float_t[fNofPtBins];
1271
1272 AliInfo("Enlarging the D0 mass windows from cut object\n");
1273 Int_t nvars = fCuts->GetNVars();
1274
1275 if(nvars<1){
1276 AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1277 return;
1278 }
1279 Float_t** rdcutsvalmine=new Float_t*[nvars];
1280 for(Int_t iv=0;iv<nvars;iv++){
1281 rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1282 }
1283
1284
1285 for (Int_t k=0;k<nvars;k++){
1286 for (Int_t j=0;j<fNofPtBins;j++){
1287
1288 // enlarge D0 window
1289 if(k==0) {
1290 fD0Window[j] =fCuts->GetCutValue(0,j);
1291 rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1292 cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1293 }
1294 else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1295
1296 // set same windows
1297 //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1298 }
1299 }
1300
1301 fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1302
1303 AliInfo("\n New windows set\n");
1304 fCuts->PrintAll();
1305
1306
1307 for(Int_t iv=0;iv<nvars;iv++){
1308 delete rdcutsvalmine[iv];
1309 }
1310 delete [] rdcutsvalmine;
1311
1312}
1313
1314
1315//____________________________ Run checks on event mixing ___________________________________________________
1316void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
5d3b1a33 1317 // check this function
5d3cf93b 1318
1319 AliCentrality *centralityObj = 0;
1320 Int_t multiplicity = -1;
1321 Double_t MultipOrCent = -1;
1322
1323 // get the pool for event mixing
1324 if(fSystem != AA){ // pp
1325 multiplicity = AOD->GetNTracks();
1326 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1327 }
1328 if(fSystem == AA){ // PbPb
1329
1330 centralityObj = AOD->GetHeader()->GetCentralityP();
1331 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1332 AliInfo(Form("Centrality is %f", MultipOrCent));
1333 }
1334
1335 AliAODVertex *vtx = AOD->GetPrimaryVertex();
1336 Double_t zvertex = vtx->GetZ(); // zvertex
1337
1338
1339
1340
1341 AliEventPool * pool = fCorrelator->GetPool();
1342
1343
1344
1345
1346 ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1347 ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1348
1349 ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1350 ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1351}
1352
1353
1354
1355
1356