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