]> 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
AliAODEvent::GetNTracks() changed with AliAODEvent::GetNumberOfTracks()
[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
c3cc0c2f 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
c3cc0c2f 261 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(7)->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
f9e355d4 420 Int_t poolbin = fAssocCuts->GetPoolBin(MultipOrCent, zVtxPosition); // get the event pool bin - commented for the moment - to be checked
5d3b1a33 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
f9e355d4 488 Int_t ncandidates = 0;
489
5d3b1a33 490 // cout << "Task debug check 3 " << endl;
491 // loop over D meson candidates
492 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
493
f9e355d4 494 //if(ncandidates) continue; // if there is more than one D candidate, skip it
495
5d3b1a33 496 // initialize variables
497 isInPeak = kFALSE;
498 isInDStarSideBand = kFALSE;
499 isInDZeroSideBand = kFALSE;
500
501 EventHasDStarCandidate = kFALSE;
502 EventHasDZeroSideBandCandidate = kFALSE;
503 EventHasDStarSideBandCandidate = kFALSE;
504
505
506 isDStarMCtag = kFALSE;
507 isDfromB = kFALSE;
508 ptDStar = -123.4;
509 phiDStar = -999;
510 etaDStar = -56.;
511 invMassDZero = - 999;
512 deltainvMDStar = -998;
513 AliAODRecoCascadeHF* dstarD0pi;
514 AliAODRecoDecayHF2Prong* theD0particle;
515 AliAODMCParticle* DStarMC=0x0;
516 Short_t daughtercharge = -2;
517 Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
518 Int_t trackiddaugh1 = -1;
519 Int_t trackidsoftPi = -1;
520 Int_t ptbin = -1;
521
522 // ============================================ using reconstruction on Data or MC
523 if(fReco){
524 // cout << "Task debug check 4 " << endl;
525 // get the D objects
526 dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
527 if(!dstarD0pi->GetSecondaryVtx()) continue;
528 theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
529 if (!theD0particle) continue;
5d3cf93b 530
5d3cf93b 531
5d3b1a33 532 // apply topological cuts
5d3cf93b 533
5d3b1a33 534 // track quality cuts
535 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
536 // region of interest + topological cuts + PID
537 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
5d3cf93b 538
5d3b1a33 539 //apply topological cuts
540 if(!isTkSelected) continue;
541 if(!isSelected) continue;
542 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
543
544 // get D candidate variables
545 ptDStar = dstarD0pi->Pt();
546 phiDStar = dstarD0pi->Phi();
547 etaDStar = dstarD0pi->Eta();
548 if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
549 if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr) efficiencyvariable = MultipOrCent;
550 if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
551 if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
552 if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
553
554
c3cc0c2f 555 // if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
5d3cf93b 556
5d3b1a33 557 phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
558 ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
5d3cf93b 559
f9e355d4 560 // cout << "DStar pt = " << ptDStar << endl;
561 // cout << "pt bin = " << ptbin << endl;
ff9e3d0b 562 // if(ptbin<3) continue;
5d3cf93b 563
5d3b1a33 564 Double_t mD0Window= fD0Window[ptbin]/3;
565 invMassDZero = dstarD0pi->InvMassD0();
566 deltainvMDStar = dstarD0pi->DeltaInvMass();
567
568
569 // get the D meson efficiency
570 DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
571
572 // check this!
573 if(fUseDmesonEfficiencyCorrection){
574 if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
c3cc0c2f 575 else DmesonWeight = 0; // do not consider te entry if the efficiency is too low
576 // else DmesonWeight = 1.;
5d3b1a33 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;
c3cc0c2f 585 // if(!isDStarMCtag) continue;
586 if(isDStarMCtag){
5d3b1a33 587 AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
588 //check if DStar from B
589 Int_t labelMother = MCDStar->GetMother();
590 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
591 if(!mother) continue;
592 Int_t motherPDG =TMath::Abs(mother->PdgCode());
593 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
c3cc0c2f 594 }
5d3b1a33 595 }
5d3cf93b 596
5d3cf93b 597
5d3b1a33 598 // fill mass histograms
ff9e3d0b 599 // cout << "crash here 1" << endl;
f9e355d4 600 // plot D0 vs DStar mass
601 if(!fmixing){
ff9e3d0b 602 cout << Form("histDZerovsDStarMass_%d",ptbin) <<endl;
f9e355d4 603 ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMass_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar);
604 if(fUseDmesonEfficiencyCorrection) ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMassWeight_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar,DmesonWeight);
605 }
ff9e3d0b 606 // cout << "crash here 2" << endl;
f9e355d4 607
608
5d3b1a33 609 // fill D0 invariant mass
610 if(!fmixing) {
5d3b1a33 611 ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
612 // cout << "Task debug check 5.1 " << endl;
613 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
614 } // end if !mixing
5d3cf93b 615
5d3b1a33 616
5d3cf93b 617
5d3cf93b 618
5d3cf93b 619
5d3b1a33 620 // good D0 canidates
621 if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
622 // fill D* invariant mass
623 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
624 // fill D* invariant mass if weighting
625 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
626 isInPeak = kTRUE;
627 // fill other good candidate distributions - pt, phi eta
628 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
629 ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
630 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
631 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
632 nOfDStarCandidates++;
f9e355d4 633 ncandidates ++;
5d3b1a33 634 EventHasDStarCandidate=kTRUE;
635 } // end if in good D* mass window
636
637 // count D* sideband candidates
638 if(fBkgMethod == kDStarSB ){
136c2bee 639 if(deltainvMDStar>0.1473 && deltainvMDStar<0.1644){
c3cc0c2f 640 // if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
5d3b1a33 641 ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
642 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
643 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
644 nOfSBCandidates++;
f9e355d4 645 ncandidates ++;
5d3b1a33 646 EventHasDStarSideBandCandidate = kTRUE;
647 }
648
649 } // end if using DStar sidebans
650
651
652 }// end good D0
653
654 // cout << "Task debug check 6 " << endl;
655 // Sideband D0
656 if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
657 isInDZeroSideBand = kTRUE;
658 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
659 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
660
661 if(fBkgMethod == kDZeroSB){
662 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
663
664 ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
665 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
666 if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
667 nOfSBCandidates++;
f9e355d4 668 ncandidates ++;
5d3b1a33 669 EventHasDZeroSideBandCandidate = kTRUE;
670 }
671 }
672 }
673 }// end SideBand D0
674 // cout << "Task debug check 7 " << endl;
5d3cf93b 675
5d3b1a33 676 if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
c3cc0c2f 677 if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.1644) continue; // skip all D* canidates outside the interesting mass ranges
5d3b1a33 678 // cout << "Good D* candidate" << endl;
5d3cf93b 679
5d3b1a33 680 // get charge of soft pion
681 daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
5d3cf93b 682
5d3b1a33 683 // get ID of daughters used for reconstruction
684 trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
685 trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
686 trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
687 }// end if reconstruction
688
ff9e3d0b 689 // cout << "crash here 3" << endl;
5d3b1a33 690
691 // ============================================ using MC kinematics only
692 Int_t DStarLabel = -1;
693
694 if(!fReco){ // use pure MC information
ff9e3d0b 695 // cout << "0 - Running on kine " << endl;
5d3b1a33 696 // get the DStar Particle
697 DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
698 if (!DStarMC) {
699 AliWarning("Careful: DStar MC Particle not found in tree, skipping");
700 continue;
701 }
ff9e3d0b 702 // cout << "1 - Have D* " << endl;
5d3b1a33 703 DStarLabel = DStarMC->GetLabel();
704 if(DStarLabel>0)isDStarMCtag = kTRUE;
5d3cf93b 705
5d3b1a33 706 Int_t PDG =TMath::Abs(DStarMC->PdgCode());
707 if(PDG !=413) continue; // skip if it is not a DStar
ff9e3d0b 708 // cout << "PDG = " << PDG << endl;
5d3b1a33 709 // check fiducial acceptance
710 if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
ff9e3d0b 711 // cout << "2 - Have D* in fiducial acceptance " << endl;
5d3cf93b 712
ff9e3d0b 713 if(DStarMC->Pt()<fCuts->GetMinPtCandidate()||DStarMC->Pt()>fCuts->GetMaxPtCandidate()) continue;
714 ptbin=fCuts->PtBin(DStarMC->Pt());
715 // cout << "3 - Have D* in ptrange acceptance " << endl;
5d3cf93b 716 //check if DStar from B
717 Int_t labelMother = DStarMC->GetMother();
718 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
5d3b1a33 719 if(!mother) continue;
5d3cf93b 720 Int_t motherPDG =TMath::Abs(mother->PdgCode());
721 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
5d3b1a33 722
723
724 Bool_t isDZero = kFALSE;
725 Bool_t isSoftPi = kFALSE;
726
727 if(fUseHadronicChannelAtKineLevel){
728 //check decay channel on MC ************************************************
5d3cf93b 729 Int_t NDaugh = DStarMC->GetNDaughters();
730 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
731
5d3b1a33 732 for(Int_t i=0; i<NDaugh;i++){ // loop on DStar daughters
733 Int_t daugh_label = DStarMC->GetDaughter(i);
734 AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
735 if(!mcDaughter) continue;
736 Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
737 if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
738
739 if(daugh_pdg == 421) {
740 Int_t NDaughD0 = mcDaughter->GetNDaughters();
741 if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
742 trackiddaugh0 = mcDaughter->GetDaughter(0);
743 trackiddaugh1 = mcDaughter->GetDaughter(1);
744 Bool_t isKaon = kFALSE;
745 Bool_t isPion = kFALSE;
746
747 for(Int_t k=0;k<NDaughD0;k++){
748 Int_t labelD0daugh = mcDaughter->GetDaughter(k);
749 AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
750 if(!mcGrandDaughter) continue;
751 Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
752 if(granddaugh_pdg==321) isKaon = kTRUE;
753 if(granddaugh_pdg==211) isPion = kTRUE;
5d3cf93b 754 }
5d3b1a33 755 if(!isKaon || !isPion) break; // skip if not correct decay channel of D0
756 isDZero = kTRUE;
757 } // end check on Dzero
758
759 if(daugh_pdg == 211) {
760 isSoftPi = kTRUE;
761 daughtercharge = mcDaughter->Charge();
762 trackidsoftPi = daugh_label;}
763 } // end loop on D* daughters
764 if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
765 } // end check decay channel
766
767 ptDStar = DStarMC->Pt();
768 phiDStar = DStarMC->Phi();
769 etaDStar = DStarMC->Eta();
770
771 if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
ff9e3d0b 772 // cout << "Dstars are selected" << endl;
5d3b1a33 773
774 }// end if pure MC information
775
776
777
778
779
5d3cf93b 780 // getting the number of triggers in the MCtag D* case
5d3b1a33 781 if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar);
782 if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
783 if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
784
785
786 fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
787 fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
788
789
ff9e3d0b 790 // cout << "crash here 4" << endl;
5d3b1a33 791
792 // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
793
f9e355d4 794 // cout << "Correlating " << endl;
5d3b1a33 795
796 Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
797
798 // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
799 if(fmixing && !execPool) {
800 AliInfo("Mixed event analysis: pool is not ready");
5d3cf93b 801 if(!isEventMixingFilledPeak && isInPeak) {
5d3b1a33 802 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
803 isEventMixingFilledPeak = kTRUE;
5d3cf93b 804 }
805 if (!isEventMixingFilledSB && isInDZeroSideBand) {
5d3b1a33 806 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
807 isEventMixingFilledSB=kTRUE;
5d3cf93b 808 }
5d3b1a33 809 continue;
810 } // end if pool not ok
811 // if analysis is on mixed event and pool settings are satisfied, fill relevant histograms and continue
812 if(fmixing&&execPool){
5d3cf93b 813 // pool is ready - run checks on bins filling
814 if(!isEventMixingFilledPeak && isInPeak) {
5d3b1a33 815 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
816 if(fFullmode) EventMixingChecks(aodEvent);
817 isEventMixingFilledPeak = kTRUE;
5d3cf93b 818 }
819
820 if(!isEventMixingFilledSB && isInDZeroSideBand) {
5d3b1a33 821 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
822 isEventMixingFilledSB=kTRUE;
5d3cf93b 823 }
5d3b1a33 824 } // end if pool ok
825
826
827
828
829 Int_t NofEventsinPool = 1;
830 if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
831
ff9e3d0b 832 Bool_t *trackOrigin = NULL;
833 // cout << "crash here 5" << endl;
5d3b1a33 834 //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
835
836 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
5d3cf93b 837
5d3b1a33 838 Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); // process the associated tracks
839 if(!analyzetracks) {
840 AliInfo("AliHFCorrelator::Cannot process the track array");
841 continue;
842 }
5d3cf93b 843
f9e355d4 844 Double_t arraytofill[5];
ff9e3d0b 845 Double_t MCarraytofill[6]; // think about this
5d3cf93b 846 Double_t weight;
847
5d3b1a33 848 Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
5d3cf93b 849
5d3b1a33 850 //************************************************** LOOP ON TRACKS *****************************************************
ff9e3d0b 851 // cout << "crash here 6" << endl;
5d3b1a33 852 for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
ff9e3d0b 853 // cout << "crash correlation 1" << endl;
5d3b1a33 854 Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
855 if(!runcorrelation) continue;
856
857 Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
858 Double_t DeltaEta = fCorrelator->GetDeltaEta();
859
860 AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
861 if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
862
c3cc0c2f 863
864 Int_t trackid = hadron->GetID();
865
866 // check D if it is a D meson daughter
867 if(!fmixing && fReco){ // for reconstruction
868 if(trackid == trackiddaugh0) continue;
869 if(trackid == trackiddaugh1) continue;
870 if(trackid == trackidsoftPi) continue;
871 }
872
873
874 Int_t label = hadron->GetLabel();
875 if(!fmixing && !fReco){ // for kinematic MC
876 AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
877 if(!part) continue;
878 if(IsDDaughter(DStarMC, part)) continue;
879 cout << "Skipping DStar daugheter " << endl;
880 }
881 // if it is ok, then do the rest
882
5d3b1a33 883 Double_t ptHad = hadron->Pt();
884 Double_t phiHad = hadron->Phi();
c3cc0c2f 885 phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
5d3b1a33 886 Double_t etaHad = hadron->Eta();
c3cc0c2f 887
888
5d3b1a33 889 Double_t efficiency = hadron->GetWeight();
c3cc0c2f 890
891
892
893
894 if(fmontecarlo) trackOrigin = fAssocCuts->IsMCpartFromHF(label,fmcArray);
ff9e3d0b 895 // cout << "crash correlation 3" << endl;
5d3b1a33 896
897 if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
898
899 ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
900 ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
901 isTrackForPeakFilled = kTRUE; // makes sure you do not fill twice in case of more candidates
902 }
903
904 if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
905 ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
906 ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
907 isTrackForSBFilled = kTRUE;
908 }
909
910 if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
911 ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
912 ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
913 isTrackForSBFilled = kTRUE;
914 }
ff9e3d0b 915 // cout << "crash correlation 4" << endl;
5d3b1a33 916
917 weight = 1;
918 if(fUseEfficiencyCorrection && efficiency){
919 weight = DmesonWeight * (1./efficiency);
920 }
921
ff9e3d0b 922 // cout << "crash correlation 5" << endl;
5d3b1a33 923 arraytofill[0] = DeltaPhi;
924 arraytofill[1] = deltainvMDStar;
925 arraytofill[2] = DeltaEta;
926 arraytofill[3] = ptHad;
f9e355d4 927 arraytofill[4] = poolbin;
5d3b1a33 928
ff9e3d0b 929 if(fmontecarlo){
930 MCarraytofill[0] = DeltaPhi;
931 MCarraytofill[1] = deltainvMDStar;
932 MCarraytofill[2] = DeltaEta;
933 MCarraytofill[3] = ptHad;
934 MCarraytofill[4] = poolbin;
935
936 if(trackOrigin[0]) MCarraytofill[5] = 1 ;
937 else if(trackOrigin[1]) MCarraytofill[5] = 2 ;
938 else MCarraytofill[5] = 0;
939 }
c3cc0c2f 940
5d3b1a33 941
942 // ============================================= FILL CORRELATION THNSparses ===============================
943
944 // filling signal
c3cc0c2f 945 // if(isInPeak){
946 if(EventHasDStarCandidate){
ff9e3d0b 947 // cout << "Filling signal " << endl;
948 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
5d3b1a33 949 //cout ("CorrelationsDStarHadron_%d",ptbin)
950 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
951 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
952 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
953 }
954
955 // filling bkg
c3cc0c2f 956 // if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
957 if(fBkgMethod == kDStarSB && EventHasDStarSideBandCandidate){ // bkg from DStar
ff9e3d0b 958 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
959 // cout << "Filling bkg from D* sidebands " << endl;
5d3b1a33 960 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
961 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
962 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
963
964 } // end if bkg from DStar
965
c3cc0c2f 966 // if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
967 if(fBkgMethod == kDZeroSB && EventHasDZeroSideBandCandidate){
5d3b1a33 968 // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
ff9e3d0b 969 // cout << "Filling bkg from Dzero sidebands " << endl;
5d3b1a33 970 if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
971 if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
972 if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
973
974 } // end if bkg from DZero
975
ff9e3d0b 976
977
978 if(fmontecarlo){
979 // add the montecarlos
980 if(!fReco && TMath::Abs(etaHad)>0.8) continue;
981 if(!isDfromB){
982 //cout << "Filling correlations from charm " << endl;
983 // cout << "Ik zoek op " << Form("CorrelationsMCfromCharmHadron_%d",ptbin) << endl;
984 // cout << "de lijst bevat : " << endl;
985 fOutputMC->ls();
986 if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
987 if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
988 if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
989 }
990 if(isDfromB){
991 //cout << "Filling correlations from beauty " << endl;
992 if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
993 if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
994 if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
995 }
996 }
997
5d3b1a33 998
999 } // end loop on associated tracks
5d3cf93b 1000
5d3b1a33 1001 } // end loop on events in event pool
5d3cf93b 1002
c3cc0c2f 1003 if(fmixing){
1004 ((TH1D*)fEMOutput->FindObject("PoolBinDistribution"))->Fill(poolbin);
1005 ((TH2D*)fEMOutput->FindObject("EventDistributionPerPoolBin"))->Fill(poolbin,NofEventsinPool);
1006 }
5d3b1a33 1007 } // end loop on D mesons
5d3cf93b 1008
f9e355d4 1009 if(!fmixing){
1010 if(nOfDStarCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfDCandidates"))->Fill(nOfDStarCandidates);
1011 if(nOfSBCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfSBCandidates"))->Fill(nOfSBCandidates);
1012 }
5d3cf93b 1013
1014
5d3b1a33 1015 Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
1016
c3cc0c2f 1017
1018
5d3b1a33 1019
1020 if(!updated) AliInfo("Pool was not updated");
1021
5d3cf93b 1022
1023} //end the exec
1024
1025//________________________________________ terminate ___________________________
1026void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
1027{
1028 // The Terminate() function is the last function to be called during
1029 // a query. It always runs on the client, it can be used to present
1030 // the results graphically or save the results to file.
1031
1032 AliAnalysisTaskSE::Terminate();
1033
1034 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1035 if (!fOutput) {
1036 printf("ERROR: fOutput not available\n");
1037 return;
1038 }
1039
1040 return;
1041}
1042//_____________________________________________________
1043Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
1044
1045 //Daughter removal in MCKine case
1046 Bool_t isDaughter = kFALSE;
1047 Int_t labelD0 = d->GetLabel();
1048
1049 Int_t mother = track->GetMother();
1050
1051 //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
1052 while (mother > 0){
1053 AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
1054 if (mcMoth){
1055 if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
1056 mother = mcMoth->GetMother(); //goes back by one
1057 } else{
1058 AliError("Failed casting the mother particle!");
1059 break;
1060 }
1061 }
1062
1063 return isDaughter;
1064}
1065
1066//_____________________________________________________
1067void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
1068
5d3cf93b 1069 Double_t Pi = TMath::Pi();
1070 Int_t nbinscorr = fPhiBins;
5d3b1a33 1071 Double_t lowcorrbin = -0.5*Pi ;
1072 Double_t upcorrbin = 1.5*Pi ;
5d3cf93b 1073
5d3cf93b 1074
5d3b1a33 1075 // create ThNSparses
5d3cf93b 1076
5d3b1a33 1077 Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
5d3cf93b 1078
1079
5d3b1a33 1080 //sparse bins
1081
1082 //1 delta_phi
1083 //2 invariant mass D *
1084 //3 delta eta
1085 //4 track pt
5d3cf93b 1086
5d3cf93b 1087
f9e355d4 1088 Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
5d3cf93b 1089
5d3cf93b 1090
5d3b1a33 1091 // for reconstruction on Data and MC
f9e355d4 1092 Int_t nbinsSparse[5]= {nbinscorr , 32 , 32, 10,nbinsPool};
1093 Double_t binLowLimitSparse[5]={lowcorrbin,0.14314 ,-1.6, 0,-0.5};
1094 Double_t binUpLimitSparse[5]= {upcorrbin ,0.14794 , 1.6, 5,nbinsPool-0.5};
5d3cf93b 1095
c3cc0c2f 1096 // Int_t nbinsSparseDStarSB[5]= {nbinscorr , 40 , 32, 10,nbinsPool};
1097 // Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.14788 ,-1.6, 0,-0.5};
1098 // Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1504 , 1.6, 5,nbinsPool-0.5};
1099
1100 Int_t nbinsSparseDStarSB[5]= {nbinscorr , 27 , 32, 10,nbinsPool};
1101 Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.1473 ,-1.6, 0,-0.5};
1102 Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1644 , 1.6, 5,nbinsPool-0.5};
5d3cf93b 1103
ff9e3d0b 1104 Int_t nbinsSparseMC[6]= {nbinscorr , 32 , 32, 10,nbinsPool,3};
1105 Double_t binLowLimitSparseMC[6]={lowcorrbin,0.14314 ,-1.6, 0,-0.5,-0.5};
1106 Double_t binUpLimitSparseMC[6]= {upcorrbin ,0.14794 , 1.6, 5,nbinsPool-0.5,2.5};
1107
5d3b1a33 1108 TString signalSparseName = "";
1109 TString bkgSparseName = "";
ff9e3d0b 1110 TString MCSparseNameCharm = "";
1111 TString MCSparseNameBeauty = "";
5d3cf93b 1112
5d3b1a33 1113 THnSparseF * CorrelationsSignal = NULL;
1114 THnSparseF * CorrelationsBkg = NULL;
5d3cf93b 1115
ff9e3d0b 1116 THnSparseF * MCCorrelationsCharm = NULL;
1117 THnSparseF * MCCorrelationsBeauty = NULL;
1118
5d3cf93b 1119
5d3b1a33 1120 Float_t * ptbinlims = fCuts->GetPtBinLimits();
1121
f9e355d4 1122
5d3cf93b 1123
5d3cf93b 1124
5d3b1a33 1125 for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1126
1127
1128
1129 if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
f9e355d4 1130
1131
5d3b1a33 1132
1133 signalSparseName = "CorrelationsDStar";
1134 if(fselect==1) signalSparseName += "Hadron_";
1135 if(fselect==2) signalSparseName += "Kaon_";
1136 if(fselect==3) signalSparseName += "KZero_";
1137
1138 bkgSparseName = "CorrelationsBkg";
1139 if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
1140 if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
1141 if(fselect==1) bkgSparseName += "Hadron_";
1142 if(fselect==2) bkgSparseName += "Kaon_";
1143 if(fselect==3) bkgSparseName += "KZero_";
1144
ff9e3d0b 1145 MCSparseNameCharm = "CorrelationsMCfromCharm";
1146 if(fselect==1) MCSparseNameCharm += "Hadron_";
1147 if(fselect==2) MCSparseNameCharm += "Kaon_";
1148 if(fselect==3) MCSparseNameCharm += "KZero_";
1149
1150 MCSparseNameBeauty = "CorrelationsMCfromBeauty";
1151 if(fselect==1) MCSparseNameBeauty += "Hadron_";
1152 if(fselect==2) MCSparseNameBeauty += "Kaon_";
1153 if(fselect==3) MCSparseNameBeauty += "KZero_";
1154
5d3b1a33 1155 signalSparseName+=Form("%d",iBin);
1156 bkgSparseName+=Form("%d",iBin);
ff9e3d0b 1157 MCSparseNameCharm+=Form("%d",iBin);
1158 MCSparseNameBeauty+=Form("%d",iBin);
5d3b1a33 1159 cout << "ThNSparses name = " << signalSparseName << endl;
1160
1161 // define thnsparses for signal candidates
f9e355d4 1162 CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1163 CorrelationsSignal->Sumw2();
5d3b1a33 1164 fCorrelationOutput->Add(CorrelationsSignal);
1165
1166 // define thnsparses for bkg candidates from DStar
1167 if(fBkgMethod == kDStarSB ){
f9e355d4 1168 CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
1169 CorrelationsBkg->Sumw2();
1170 fCorrelationOutput->Add(CorrelationsBkg);
5d3b1a33 1171 }
1172
1173 // define thnsparses for bkg candidates from DZero
1174 if(fBkgMethod == kDZeroSB ){
f9e355d4 1175 CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
5d3b1a33 1176 CorrelationsBkg->Sumw2();
1177 fCorrelationOutput->Add(CorrelationsBkg);
1178 }
1179
ff9e3d0b 1180 if(fmontecarlo){
1181 MCCorrelationsCharm = new THnSparseF(MCSparseNameCharm.Data(),"Correlations for DStar from charm; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
1182 MCCorrelationsCharm->Sumw2();
1183 fOutputMC->Add(MCCorrelationsCharm);
1184
1185 MCCorrelationsBeauty = new THnSparseF(MCSparseNameBeauty.Data(),"Correlations for DStar from beauty; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
1186 MCCorrelationsBeauty->Sumw2();
1187 fOutputMC->Add(MCCorrelationsBeauty);
1188 }
1189
5d3b1a33 1190 } // end loop on bins
5d3cf93b 1191
5d3cf93b 1192
1193
1194}
5d3b1a33 1195
5d3cf93b 1196//__________________________________________________________________________________________________
1197void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
5d3b1a33 1198
1199 Double_t Pi = TMath::Pi();
5d3cf93b 1200 Int_t nbinscorr = fPhiBins;
5d3b1a33 1201 Double_t lowcorrbin = -0.5*Pi ;
1202 Double_t upcorrbin = 1.5*Pi ;
1203
1204 // event counter
1205 TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
5d3cf93b 1206 NofEvents->GetXaxis()->SetBinLabel(1," All events");
1207 NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1208 NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1209 NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1210 NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1211 NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1212 NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1213 NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1214 NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1215 NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1216 NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1217 fOutput->Add(NofEvents);
5d3cf93b 1218
5d3b1a33 1219 //event properties
1220 TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1221 fOutput->Add(EventPropsCheck);
5d3cf93b 1222
5d3b1a33 1223 //event properties
1224 TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",1000,0,1000);
1225 fOutput->Add(SPDMultiplicty);
5d3cf93b 1226
5d3cf93b 1227
5d3b1a33 1228 // =================================================== D* histograms
1229 TH1F * D0mass = NULL;
1230 TH1F * DStarMass = NULL;
1231 TH1F * DStarFromSBMass = NULL;
f9e355d4 1232 TH2D * DZerovsDStarMass = NULL;
5d3cf93b 1233
5d3b1a33 1234 TH1F * D0massWeighted = NULL;
1235 TH1F * DStarMassWeighted = NULL;
1236 TH1F * DStarFromSBMassWeighted = NULL;
f9e355d4 1237 TH2D * DZerovsDStarMassWeighted = NULL;
1238
5d3cf93b 1239
f9e355d4 1240 TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "", nameDZerovsDStarMass = "";
5d3cf93b 1241
5d3b1a33 1242 Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1243 Float_t * ptbinlims = fCuts->GetPtBinLimits();
5d3cf93b 1244
5d3b1a33 1245 //GetMinPtCandidate()
f9e355d4 1246
5d3cf93b 1247
5d3b1a33 1248 for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1249
f9e355d4 1250
5d3b1a33 1251
1252 if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1253
1254
1255 std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
1256
1257 nameDZeroMass = "histDZeroMass_";
1258 nameDStarMass = "histDStarMass_";
1259 nameDStarFromSBMass = "histDStarFromSBMass_";
f9e355d4 1260 nameDZerovsDStarMass = "histDZerovsDStarMass_";
5d3b1a33 1261
1262 nameDZeroMass+=Form("%d",iBin);
1263 nameDStarMass+=Form("%d",iBin);
1264 nameDStarFromSBMass+=Form("%d",iBin);
f9e355d4 1265 nameDZerovsDStarMass+=Form("%d",iBin);
1266 cout << "D vs D histogram: " << nameDZerovsDStarMass << endl;
5d3b1a33 1267
f9e355d4 1268 D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1269 DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1270 DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1271 DZerovsDStarMass = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
5d3b1a33 1272
1273 if(!fmixing){
f9e355d4 1274 fDmesonOutput->Add(D0mass);
1275 fDmesonOutput->Add(DStarMass);
1276 fDmesonOutput->Add(DStarFromSBMass);
1277 fDmesonOutput->Add(DZerovsDStarMass);
5d3b1a33 1278 }
1279
1280 // if using D meson efficiency, define weighted histos
1281 if(fUseDmesonEfficiencyCorrection){
f9e355d4 1282
5d3b1a33 1283 nameDZeroMass = "histDZeroMassWeight_";
1284 nameDStarMass = "histDStarMassWeight_";
1285 nameDStarFromSBMass = "histDStarFromSBMassWeight_";
f9e355d4 1286 nameDZerovsDStarMass = "histDZerovsDStarMassWeight_";
5d3b1a33 1287
1288 nameDZeroMass+=Form("%d",iBin);
1289 nameDStarMass+=Form("%d",iBin);
1290 nameDStarFromSBMass+=Form("%d",iBin);
f9e355d4 1291 nameDZerovsDStarMass+=Form("%d",iBin);
5d3b1a33 1292
f9e355d4 1293 D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1294 DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1295 DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1296 DZerovsDStarMassWeighted = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
5d3b1a33 1297
1298 if(!fmixing){
f9e355d4 1299 fDmesonOutput->Add(D0massWeighted);
1300 fDmesonOutput->Add(DStarMassWeighted);
1301 fDmesonOutput->Add(DStarFromSBMassWeighted);
1302 fDmesonOutput->Add(DZerovsDStarMassWeighted);
5d3b1a33 1303 }
1304 }
1305 }// end loop on pt bins
5d3cf93b 1306
5d3cf93b 1307
5d3b1a33 1308 // pt distributions
1309 TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
1310 fDmesonOutput->Add(RecoPtDStar);
1311 TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
1312 fDmesonOutput->Add(RecoPtBkg);
5d3cf93b 1313
f9e355d4 1314 TH1D *NumberOfDCandidates = new TH1D("NumberOfDCandidates","Number of D candidates",10,-0.5,9.5);
1315 TH1D *NumberOfSBCandidates = new TH1D("NumberOfSBCandidates","Number of SB candidates",10,-0.5,9.5);
1316 if(!fmixing) fDmesonOutput->Add(NumberOfDCandidates);
1317 if(!fmixing) fDmesonOutput->Add(NumberOfSBCandidates);
1318
5d3b1a33 1319 // phi distribution
1320 TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1321 TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
5d3cf93b 1322
5d3b1a33 1323 // eta distribution
1324 TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1325 TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
5d3cf93b 1326
5d3b1a33 1327 if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
1328 if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
1329 if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
1330 if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
5d3cf93b 1331
5d3cf93b 1332
5d3b1a33 1333 // single track related histograms
1334 // phi distribution
1335 TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1336 TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
5d3cf93b 1337
5d3b1a33 1338 // eta distribution
1339 TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1340 TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
5d3cf93b 1341
f9e355d4 1342 TH1D * TracksPerDcandidate = new TH1D("TracksPerDcandidate","Distribution of number of tracks per D meson candidate; N tracks; counts",200,-0.5,199.5);
1343 TH1D * TracksPerSBcandidate = new TH1D("TracksPerSBcandidate","Distribution of number of tracks per sideband candidate; N tracks; counts",200,-0.5,199.5);
1344
5d3b1a33 1345 if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
1346 if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
1347 if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
1348 if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
f9e355d4 1349 if(!fmixing) fTracksOutput->Add(TracksPerDcandidate);
1350 if(!fmixing) fTracksOutput->Add(TracksPerSBcandidate);
1351
5d3cf93b 1352
5d3b1a33 1353 // Montecarlo for D*
1354 TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1355 if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
5d3cf93b 1356
5d3b1a33 1357 TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1358 if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
5d3cf93b 1359
5d3b1a33 1360 TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1361 if(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
5d3cf93b 1362
5d3cf93b 1363
f9e355d4 1364
5d3b1a33 1365 // event mixing histograms
5d3cf93b 1366 TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1367 CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1368 CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1369 CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1370 CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1371
5d3b1a33 1372 if(fmixing) fEMOutput->Add(CheckPoolReadiness);
5d3cf93b 1373
5d3b1a33 1374
c3cc0c2f 1375 Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
f9e355d4 1376 Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
1377 Int_t nPoolBins = NofCentBins*NofZVrtxBins;
c3cc0c2f 1378
1379 Int_t maxevents = fAssocCuts->GetMaxNEventsInPool();
f9e355d4 1380
1381
c3cc0c2f 1382 TH1D * PoolBinDistribution = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5);
1383 if(fmixing) fEMOutput->Add(PoolBinDistribution);
f9e355d4 1384
c3cc0c2f 1385 TH2D * EventDistributionPerPoolBin = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5,maxevents+2,0,maxevents+2);
1386 if(fmixing) fEMOutput->Add(EventDistributionPerPoolBin);
1387
5d3cf93b 1388}
1389
5d3b1a33 1390
5d3cf93b 1391//__________________________________________________________________________________________________
1392void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1393
1394
1395 //Float_t* ptbins = fCuts->GetPtBinLimits();
1396 if(fD0Window) delete fD0Window;
1397 fD0Window = new Float_t[fNofPtBins];
1398
1399 AliInfo("Enlarging the D0 mass windows from cut object\n");
1400 Int_t nvars = fCuts->GetNVars();
1401
1402 if(nvars<1){
1403 AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1404 return;
1405 }
1406 Float_t** rdcutsvalmine=new Float_t*[nvars];
1407 for(Int_t iv=0;iv<nvars;iv++){
1408 rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1409 }
1410
1411
1412 for (Int_t k=0;k<nvars;k++){
1413 for (Int_t j=0;j<fNofPtBins;j++){
1414
1415 // enlarge D0 window
1416 if(k==0) {
1417 fD0Window[j] =fCuts->GetCutValue(0,j);
1418 rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1419 cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1420 }
1421 else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1422
1423 // set same windows
1424 //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1425 }
1426 }
1427
1428 fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1429
1430 AliInfo("\n New windows set\n");
1431 fCuts->PrintAll();
1432
1433
1434 for(Int_t iv=0;iv<nvars;iv++){
1435 delete rdcutsvalmine[iv];
1436 }
1437 delete [] rdcutsvalmine;
1438
1439}
1440
1441
1442//____________________________ Run checks on event mixing ___________________________________________________
1443void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
5d3b1a33 1444 // check this function
5d3cf93b 1445
1446 AliCentrality *centralityObj = 0;
1447 Int_t multiplicity = -1;
1448 Double_t MultipOrCent = -1;
1449
1450 // get the pool for event mixing
1451 if(fSystem != AA){ // pp
4640c275 1452 multiplicity = AOD->GetNumberOfTracks();
5d3cf93b 1453 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1454 }
1455 if(fSystem == AA){ // PbPb
1456
1457 centralityObj = AOD->GetHeader()->GetCentralityP();
1458 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1459 AliInfo(Form("Centrality is %f", MultipOrCent));
1460 }
1461
1462 AliAODVertex *vtx = AOD->GetPrimaryVertex();
1463 Double_t zvertex = vtx->GetZ(); // zvertex
1464
1465
1466
1467
1468 AliEventPool * pool = fCorrelator->GetPool();
1469
1470
1471
1472
1473 ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1474 ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1475
1476 ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1477 ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1478}
1479
1480
1481
1482
1483