Removed debug printout + fix in the AddTask for D* correlations (Sandro)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskDStarCorrelations.cxx
CommitLineData
815e6bab 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
28/* $Id$ */
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 "AliEventPoolManager.h"
44#include "AliVParticle.h"
45#include "AliAnalysisManager.h"
46#include "AliAODInputHandler.h"
47#include "AliAODHandler.h"
48#include "AliESDtrack.h"
49#include "AliAODMCParticle.h"
50#include "AliNormalizationCounter.h"
51#include "AliReducedParticle.h"
52
53
54
55ClassImp(AliAnalysisTaskDStarCorrelations)
56
57
58//__________________________________________________________________________
59AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
60AliAnalysisTaskSE(),
61fhandler(0x0),
62fPoolMgr(0x0),
63fmcArray(0x0),
64fCounter(0x0),
65fselect(0),
66fmontecarlo(kFALSE),
67fmixing(kFALSE),
68fEvents(0),
69fDebug(0),
70
71fOutput(0x0),
72fCuts(0),
73fCuts2(0)
74{
75// default constructor
76}
77
78//__________________________________________________________________________
79AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts) :
80AliAnalysisTaskSE(name),
81
82fhandler(0x0),
83fPoolMgr(0x0),
84fmcArray(0x0),
85fCounter(0x0),
86fselect(0),
87fmontecarlo(kFALSE),
88fmixing(kFALSE),
89fEvents(0),
90fDebug(0),
91
92fOutput(0x0),
93fCuts(0),
94fCuts2(AsscCuts)
95{
96 fCuts=cuts;
97 Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
98 DefineInput(0, TChain::Class());
99 DefineOutput(1,TList::Class()); // histos from data
100 DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
101 DefineOutput(3,AliNormalizationCounter::Class()); // normalization
102 DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my cuts
103}
104
105//__________________________________________________________________________
106
107AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
108 //
109 // destructor
110 //
111
112 Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");
113
114 if(fhandler) {delete fhandler; fhandler = 0;}
115 if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}
116 if(fmcArray) {delete fmcArray; fmcArray = 0;}
117 if(fCounter) {delete fCounter; fCounter = 0;}
118 if(fOutput) {delete fOutput; fOutput = 0;}
119 if(fCuts) {delete fCuts; fCuts = 0;}
120 if(fCuts2) {delete fCuts2; fCuts2=0;}
121
122}
123
124//___________________________________________________________
125void AliAnalysisTaskDStarCorrelations::Init(){
126 //
127 // Initialization
128 //
129 if(fDebug > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
130
131 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
132 // Post the D* cuts
133 PostData(2,copyfCuts);
134
135 // Post the hadron cuts
136 PostData(4,fCuts2);
137
138 return;
139}
140
141
142//_________________________________________________
143void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
144 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
145
146 //slot #1
147 //OpenFile(0);
148 fOutput = new TList();
149 fOutput->SetOwner();
150
151 // define histograms
152 DefineHistoForAnalysis();
153
154 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(3)->GetContainer()->GetName()));
155 fCounter->Init();
156
157
158 // definition of the Pool Manager for Event Mixing
159
160 Int_t MaxNofEvents = 200;
161 Int_t MinNofTracks = 1000;
162
163 Int_t NofMultiplicityBins = 5;
164 Double_t MBins[]={0,20,40,60,80,500};
165 Double_t * MultiplicityBins = MBins;
166
167 Int_t NofZVrtxBins = 5;
168 Double_t ZBins[]={-10,-5,-2.5,2.5,5,10};
169 Double_t *ZVrtxBins = ZBins;
170
171
172 fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofMultiplicityBins, MultiplicityBins, NofZVrtxBins, ZVrtxBins);
173
174
175}
176//_________________________________________________
177void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
178
179 cout << " " << endl;
180 cout << "=================================================================================" << endl;
181
bf658d5c 182 if(fselect==1) cout << "TASK::Correlation with hadrons "<< endl;
183 if(fselect==2) cout << "TASK::Correlation with kaons "<< endl;
184 if(fselect==3) cout << "TASK::Correlation with kzeros "<< endl;
815e6bab 185
186 if (!fInputEvent) {
187 Error("UserExec","NO EVENT FOUND!");
188 return;
189 }
190
bf658d5c 191 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
192 if(!aodEvent){
193 AliError("AOD event not found!");
194 return;
195 }
815e6bab 196 Double_t pi = TMath::Pi();
197
198 fEvents++; // event counter
199 ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
815e6bab 200 fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
bf658d5c 201 if(fmontecarlo && !fmcArray){
202 AliError("Array of MC particles not found");
203 return;
204 }
205
815e6bab 206 // initialize the pool for event mixing
207 Int_t multiplicity = aodEvent->GetNTracks();
208 AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
209 Double_t zvertex = vtx->GetZ();
210 Double_t multip = multiplicity;
211
212 if(TMath::Abs(zvertex)>=10 || multip>500 || multip == 0) {
213 AliInfo(Form("Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,multip));
214 return;
215 }
216
217
218 AliEventPool* pool = fPoolMgr->GetEventPool(multip, zvertex);
219 if (!pool) AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f cm", multip, zvertex));
220
221
222
223 // D* reconstruction
224
225 TClonesArray *arrayDStartoD0pi=0;
226
227
228 if(!aodEvent && AODEvent() && IsStandardAOD()) {
229 // In case there is an AOD handler writing a standard AOD, use the AOD
230 // event in memory rather than the input (ESD) event.
231 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
232 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
233 // have to taken from the AOD event hold by the AliAODExtension
234 AliAODHandler* aodHandler = (AliAODHandler*)
235 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
236 if(aodHandler->GetExtensions()) {
237 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
238 AliAODEvent *aodFromExt = ext->GetAOD();
239 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
240 }
241 } else {
242 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
243 }
244
245 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
246
247 // initialize variables you will need for the D*
248
249 Double_t ptDStar;//
250 Double_t phiDStar;//
251 Double_t etaDStar;//
252 Bool_t isInPeak, isInSideBand, isDStarMCtag;
253 Double_t invMassDZero;
254 Double_t deltainvMDStar;
255
256
257 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
258 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
259
260
261
262 if(fselect ==3){// check the K0 invariant mass
263 TObjArray * fillkzerospectra = AcceptAndReduceKZero(aodEvent, 0,0);
264 delete fillkzerospectra;
265 }
266
267 //MC tagging for DStar
268 //D* and D0 prongs needed to MatchToMC method
269 Int_t pdgDgDStartoD0pi[2]={421,211};
270 Int_t pdgDgD0toKpi[2]={321,211};
271
272 //loop on D* candidates
273 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
274 isInPeak = kFALSE;
275 isInSideBand = kFALSE;
276 isDStarMCtag = kFALSE;
277 ptDStar = -123.4;
278 phiDStar = -999;
279 etaDStar = -56.;
280 invMassDZero = - 999;
281 deltainvMDStar = -998;
282
283
284
285 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
286 if(!dstarD0pi->GetSecondaryVtx()) continue;
287 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
288 if (!theD0particle) continue;
289
290 // track quality cuts
291 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
292 // region of interest + topological cuts + PID
293 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
294
295 //apply selections
296 if(!isTkSelected) continue;
297 if(!isSelected) continue;
298 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
299
300 if(fmontecarlo){
301 // find associated MC particle for D* ->D0toKpi
302 Int_t mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
303 if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
304 }
305
306 ptDStar = dstarD0pi->Pt();
307 phiDStar = dstarD0pi->Phi();
308 etaDStar = dstarD0pi->Eta();
309
310 Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
311
312 Double_t dmDStarWindow =0.0019;// 0.0019 = 3 sigma
313 Double_t mD0Window=0.074;
314
315 if (ptbin==1) mD0Window = 0.026; //0.5-1
316 if (ptbin==2) mD0Window = 0.022; //1-2
317 if (ptbin==3) mD0Window = 0.024; //2-3
318 if (ptbin==4) mD0Window = 0.032;
319 if (ptbin==5) mD0Window = 0.032;
320 if (ptbin==6) mD0Window = 0.036;
321 if (ptbin==7) mD0Window = 0.036;
322 if (ptbin==8) mD0Window = 0.036;
323 if (ptbin==9) mD0Window = 0.058;
324 if (ptbin==10) mD0Window = 0.058;
325 if (ptbin>10) mD0Window = 0.074;
326
327
328 invMassDZero = dstarD0pi->InvMassD0();
329 ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
330
331 deltainvMDStar = dstarD0pi->DeltaInvMass();
332
333
334
335 //good candidates
336 if (TMath::Abs(invMassDZero-mPDGD0)<mD0Window){
337
338 ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
339 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
340
341 ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar);
342 isInPeak = kTRUE;
343 }
344 }// end if good candidates
345
346 //sidebands
347 if (TMath::Abs(invMassDZero-mPDGD0)>1.3*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<4.*mD0Window ){
348 ((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
349 ((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero);
350
351 if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
352 ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar);
353 isInSideBand = kTRUE;
354 }
355
356 }//end if sidebands
357
358
359
360 if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
361
362
363 Int_t trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
364 Int_t trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
365 Int_t trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
366
367 ptDStar = dstarD0pi->Pt();
368 phiDStar = dstarD0pi->Phi();
369 etaDStar = dstarD0pi->Eta();
370
371 if(!fmixing){ // analysis on Single Event
372
373 TObjArray* selectedtracks = new TObjArray();
374 if(fselect==1 || fselect ==2) selectedtracks = AcceptAndReduceTracks(aodEvent);
2faf0b4b 375 if(fselect==3) {selectedtracks = AcceptAndReduceKZero(aodEvent,iDStartoD0pi,1);}
815e6bab 376 Int_t noftracks = selectedtracks->GetEntriesFast();
377 Int_t counterPeak =0;
378 Int_t counterSB = 0;
379
380 for(Int_t i =0; i<noftracks; i++){ // loop on tracks/k0 candidates in aod event
381
382 AliReducedParticle *redpart = (AliReducedParticle*)selectedtracks->At(i);
383 Double_t phiHad=redpart->Phi();
384
385 Double_t ptHad=redpart->Pt();
386 Double_t etaHad=redpart->Eta();
387
388 Int_t label = redpart->GetLabel();
389
390 Int_t trackid = redpart->GetID();
391
392
393 // check that the track is not a D* daughter
394 if(trackid == trackiddaugh0) continue;
395 if(trackid == trackiddaugh1) continue;
396 if(trackid == trackidsoftPi) continue;
397
398 if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
399
400 Int_t PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
401
402 FillMCTagCorrelations(ptDStar,phiDStar,etaDStar,ptHad,phiHad,etaHad,PartSource);
403 }
404 if(isInPeak) {
405 FillCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad);// function for correlations
406 counterPeak++;
407 if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
408 if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
409
410 ((TH1F*)fOutput->FindObject("PhiTrigger"))->Fill(phiDStar);
411
412
413
414 if (phiHad > 1.5*pi) phiHad = phiHad - 2*pi;
415 if (phiHad < -0.5*pi) phiHad = phiHad + 2*pi;
416 ((TH1F*)fOutput->FindObject("PhiPart"))->Fill(phiHad);
417
418 }
2faf0b4b 419
815e6bab 420 if(isInSideBand) {
421
422 FillSideBandCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad); // function for sidebands
423 if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
424 if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
425 ((TH1F*)fOutput->FindObject("PhiSideBand"))->Fill(phiDStar);
426
427 counterSB++;
428 }
2faf0b4b 429
815e6bab 430 } // end loop on tracks
431
432 if(counterPeak) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(counterPeak);
433 if(counterSB) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(counterSB);
434
435
436
437
438 } // end if SE Analysis
439
440 if(fmixing) { // analysis on Mixed Events
441 if (pool->IsReady()|| pool->GetCurrentNEvents()>=5){ // check if the pool is ready
442
443 pool->PrintInfo();
444
445 Int_t multbinflag = pool->MultBinIndex();
446 Int_t zvtxflag = pool->ZvtxBinIndex();
447
2faf0b4b 448 if(isInPeak) ((TH2I*)fOutput->FindObject("EventMixingCheck"))->Fill(multbinflag,zvtxflag);
449
815e6bab 450
bf658d5c 451 TObjArray* mixedtracks = 0x0;
815e6bab 452
453 for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {//loop over the events in the pool
454 mixedtracks = pool->GetEvent(jMix);
455 if(!mixedtracks) cout << "No Mixed tracks " << endl;
456 Int_t jMax = mixedtracks->GetEntriesFast();
457
458 for(Int_t iMix =0; iMix<jMax; iMix++){ //loop on the tracks of the event
459 AliVParticle *redpart = (AliVParticle*)mixedtracks->At(iMix);
460 Double_t phiHad=redpart->Phi();
461 Double_t etaHad=redpart->Eta();
462 Double_t ptHad=redpart->Pt();
463 Int_t label = redpart->GetLabel();
464
465
466 if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
467
468 Int_t PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
469
470 FillMCTagCorrelations(ptDStar,phiDStar,etaDStar,ptHad,phiHad,etaHad,PartSource);
471 }
472
473 if(isInPeak) {
474 FillCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad);// function for correlations
475
476 if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
477 if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
478
479 ((TH1F*)fOutput->FindObject("PhiTrigger"))->Fill(phiDStar);
480 if (phiHad > 1.5*pi) phiHad = phiHad - 2*pi;
481 if (phiHad < -0.5*pi) phiHad = phiHad + 2*pi;
482 ((TH1F*)fOutput->FindObject("PhiPart"))->Fill(phiHad);
483
484
485
486 }
487
488 if(isInSideBand) FillSideBandCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad); // function for sidebands
489
490 } // end loop on tracks
491 }// end loop on events in pool
492 } // end if pool is ready
493
494 } // end ME analysis
495
496 }// end loop on D* candidates
497
498 if(fmixing) { // update the pool for Event Mixing
499 if(fselect==1 || fselect==2)pool->UpdatePool(AcceptAndReduceTracks(aodEvent)); // updating the pool for hadrons
500 if(fselect==3) pool->UpdatePool(AcceptAndReduceKZero(aodEvent,0,0)); // updating the pool for K0s
501 }
502 //cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END OF THE EVENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
503
504 PostData(1,fOutput); // set the outputs
505 PostData(3,fCounter); // set the outputs
506
507} //end the exec
508
509//________________________________________ terminate ___________________________
510void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
511{
512 // The Terminate() function is the last function to be called during
513 // a query. It always runs on the client, it can be used to present
514 // the results graphically or save the results to file.
515
516 AliAnalysisTaskSE::Terminate();
517
518 fOutput = dynamic_cast<TList*> (GetOutputData(1));
519 if (!fOutput) {
520 printf("ERROR: fOutput not available\n");
521 return;
522 }
523
524 return;
525}
526
527
528//_____________________________________________________
529TObjArray* AliAnalysisTaskDStarCorrelations::AcceptAndReduceTracks(AliAODEvent* inputEvent){
530
531 Int_t nTracks = inputEvent->GetNTracks();
532 AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
533 Double_t Bz = inputEvent->GetMagneticField();
534
535
536 TObjArray* tracksClone = new TObjArray;
537 tracksClone->SetOwner(kTRUE);
538 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
539 AliAODTrack* track = inputEvent->GetTrack(iTrack);
540 if (! track) continue;
541
542 if(!fCuts2->IsHadronSelected(track,vtx,Bz)) continue; // apply selection of hadrons
543
544
545 if(fselect ==2){
546 if(!fCuts2->CheckKaonCompatibility(track,fmontecarlo,fmcArray)) continue; // check if it is a Kaon - data and MC
547 }
548
549 AliVParticle * part = (AliVParticle*)track;
550 tracksClone->Add(new AliReducedParticle(part->Eta(), part->Phi(), part->Pt(),track->GetLabel(),track->GetID()));
551
552 }
553 return tracksClone;
554}
555
556//_____________________________________________________
557TObjArray* AliAnalysisTaskDStarCorrelations::AcceptAndReduceKZero(AliAODEvent* inputEvent, Int_t loopindex, Int_t plotassociation){
558
559 Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
560 TObjArray* KZeroClone = new TObjArray;
561 AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();
562 Int_t v0label = -1;
563 Int_t pdgDgK0toPipi[2] = {211,211};
564 Double_t mPDGK0=TDatabasePDG::Instance()->GetParticle(310)->Mass();
565 const Int_t size = inputEvent->GetNumberOfV0s();
566 Int_t prevnegID[size];
567 Int_t prevposID[size];
568 for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
569 AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);
570 if(!v0) {cout << "is not a v0 at step " << iv0 << endl; continue;}
571
572 if(!fCuts2->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
573
574 // checks to avoid double counting
575 Int_t negID = -999;
576 Int_t posID = -998;
577 //Int_t a = 0;
578 prevnegID[iv0] = -997;
579 prevposID[iv0]= -996;
580 negID = v0->GetNegID();
581 posID = v0->GetPosID();
582 Bool_t DoubleCounting = kFALSE;
583
584 for(Int_t k=0; k<iv0; k++){
585 if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){
586 DoubleCounting = kTRUE;
587 //a=k;
588 break;
589 }//end if
590 } // end for
591
592 if(DoubleCounting) {cout << "SKIPPING DOUBLE COUNTING" << endl;continue;}
593 else{ prevposID[iv0] = posID; prevnegID[iv0] = negID;}
594
595 if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
596 Double_t k0pt = v0->Pt();
597 Double_t k0eta = v0->Eta();
598 Double_t k0Phi = v0->Phi();
599 Double_t k0InvMass = v0->MassK0Short();
600
601 if (loopindex == 0) {
602 if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
603 if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*
604 }
605 // if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
606
607 if(TMath::Abs(k0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
608 KZeroClone->Add(new AliReducedParticle(k0eta , k0Phi, k0pt,v0label,0));
609
610 }
611
612 return KZeroClone;
613}
614
615//_____________________________________________________
616void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
617
618 Double_t Pi = TMath::Pi();
619 Int_t nbinscorr = 32;
620 Double_t lowcorrbin = -0.5*Pi - Pi/32; // shift the bin by half the width so that at 0 is it the bin center
621 Double_t upcorrbin = 1.5*Pi - Pi/32;
622
623 // ========================= histograms for both Data and MonteCarlo
624
625
626 TH1D * NofEvents = new TH1D("NofEvents","NofEvents",1,0,1);
627 fOutput->Add(NofEvents);
628
629 TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
630 fOutput->Add(D0InvMass);
631
632 TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);
633 fOutput->Add(D0InvMassinSB);
634
635 TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
636 fOutput->Add(DeltaInvMass);
637
638 TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
639 fOutput->Add(bkgDeltaInvMass);
640
641 TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",30,0,30);
642 fOutput->Add(RecoPtDStar);
643
644 TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",30,0,30);
645 fOutput->Add(RecoPtBkg);
646
647 TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
648 if(fselect==3) fOutput->Add(KZeroSpectra);
649
650 TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);
651 if(fselect==3) fOutput->Add(KZeroSpectraifHF);
652
653 TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","NofTracksInPeak",500,0.5,500.5);
654 fOutput->Add(NofTracksInPeak);
655
656 TH1D * NofTracksInSB = new TH1D("NofTracksInSB","NofTracksInSB",500,0.5,500.5);
657 fOutput->Add(NofTracksInSB);
658
659 TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
660 if(fmixing) fOutput->Add(EventMixingCheck);
661
662
663
664 TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
665 MCSources->GetXaxis()->SetBinLabel(1,"All ");
666 MCSources->GetXaxis()->SetBinLabel(2," from Heavy flavour");
667 MCSources->GetXaxis()->SetBinLabel(3," from c->D");
668 MCSources->GetXaxis()->SetBinLabel(4," from b->D");
669 MCSources->GetXaxis()->SetBinLabel(5," from b->B");
670 if(fmontecarlo) fOutput->Add(MCSources);
671
672 TH1F * PhiTrigger = new TH1F("PhiTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi);
673 fOutput->Add(PhiTrigger);
674
675 TH1F * PhiSideBand = new TH1F("PhiSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi);
676 fOutput->Add(PhiSideBand);
677
678 TH1F * PhiPart = new TH1F("PhiPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi);
679 fOutput->Add(PhiPart);
680
681
682 //correlations histograms
683 TString histoname1 = "DPhiDStar";
684 if(fselect==1) histoname1 += "Hadron";
685 if(fselect==2) histoname1 += "Kaon";
686 if(fselect==3) histoname1 += "KZero";
687
688
689 TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
690
691 TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
692
693 //side band background histograms
694 TString histoname2 = "bkg";
695 histoname2 += histoname1;
696 TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
697 TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
698
699
700 fOutput->Add(DPhiDStar);
701
702 if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
703
704 fOutput->Add(bkgDPhiDStar);
705
706 if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
707
708
709 // ========================= histos for analysis on MC
710 TString histoname3 = "MCTag";
711 histoname3 += histoname1;
712 TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
713
714 TString histoname44 = "CharmDOrigin";
715 histoname44 += histoname1;
716 histoname44 += "MC";
717
718 TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
719
720
721 TString histoname54 = "BeautyDOrigin";
722 histoname54 += histoname1;
723 histoname54 += "MC";
724 TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
725
726 TString histoname55 = "BeautyBOrigin";
727 histoname55 += histoname1;
728 histoname55 += "MC";
729 TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
730
731 if(fmontecarlo){
732
733 fOutput->Add(MCTagDPhiDStar);
734 fOutput->Add(CharmDOriginDPhiDStar);
735 fOutput->Add(BeautyDOriginDPhiDStar);
736 fOutput->Add(BeautyBOriginDPhiDStar);
737
738 }
739
740
741
742}
743
744
745//____________________________ Function for correlations ___________________________________________________
746void AliAnalysisTaskDStarCorrelations::FillCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack){
747 Double_t pi = TMath::Pi();
748 Double_t deltaPhi, deltaEta;
749 deltaPhi = phiTrig - phiTrack;
750 deltaEta = etaTrig - etaTrack;
751 // set correct Delta Phi range
752 if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
753 if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
2faf0b4b 754
815e6bab 755 if(fselect==1) ((TH3D*)fOutput->FindObject("DPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
756 if(fselect==2) ((TH3D*)fOutput->FindObject("DPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
757 if(fselect==3) ((TH3D*)fOutput->FindObject("DPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
2faf0b4b 758
815e6bab 759
760 return;
761}
762
763//____________________________ Function for sidebands ___________________________________________________
764void AliAnalysisTaskDStarCorrelations::FillSideBandCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack){
765
766 Double_t pi = TMath::Pi();
767 Double_t deltaPhi, deltaEta;
768 deltaPhi = phiTrig - phiTrack;
769 deltaEta = etaTrig - etaTrack;
770 // set correct Delta Phi range
771 if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
772 if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
2faf0b4b 773
815e6bab 774 if(fselect==1) ((TH3D*)fOutput->FindObject("bkgDPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
775 if(fselect==2) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
776 if(fselect==3) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
2faf0b4b 777
815e6bab 778
779 return;
780
781
782}
783
784//____________________________ Function for sidebands ___________________________________________________
785void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t ptTrack, Double_t phiTrack, Double_t etaTrack, Int_t mcSource){
786Double_t deltaPhi = phiTrig - phiTrack;
787Double_t deltaEta = etaTrig - etaTrack;
788// set correct Delta Phi range
789
790 Double_t pi = TMath::Pi();
791
792if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
793if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
2faf0b4b 794
815e6bab 795 if(fselect==1) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
796 if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
797 if(fselect==3) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
798
799
800
801((TH1F*)fOutput->FindObject("MCSources"))->Fill(0);
2faf0b4b 802
815e6bab 803if (mcSource==44){ // is from charm ->D
804 if(fselect==1) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
805 if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
806 if(fselect==3) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
807
808
809 ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
810 ((TH1F*)fOutput->FindObject("MCSources"))->Fill(2);
811 }
2faf0b4b 812
815e6bab 813if (mcSource==54){ // is from beauty -> D
814 if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
815 if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
816 if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
817 if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
818 if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(3);
819 }
2faf0b4b 820
815e6bab 821if (mcSource==55){ // is from beauty ->B
822 if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
823 if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
824 if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyOriginBDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
825 if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
826 if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(4);
827 }
828 return;
829}
830
831
832
833
834