Fixes for std:: need with the trunk of Root (Jochen, Yves)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliHFCorrelator.cxx
CommitLineData
bce70c96 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 Heavy Flavour Correlations Analysis
18// Single Event and Mixed Event Analysis are implemented
19//
20//-----------------------------------------------------------------------
21//
22//
23// Author S.Bjelogrlic
24// Utrecht University
25// sandro.bjelogrlic@cern.ch
26//
27//-----------------------------------------------------------------------
28
29/* $Id$ */
30
31#include <TParticle.h>
32#include <TVector3.h>
33#include <TChain.h>
34#include "TROOT.h"
35#include "AliHFCorrelator.h"
36#include "AliRDHFCutsDStartoKpipi.h"
37#include "AliHFAssociatedTrackCuts.h"
38#include "AliEventPoolManager.h"
39#include "AliReducedParticle.h"
40#include "AliCentrality.h"
41
e61afb80 42using std::cout;
43using std::endl;
44
bce70c96 45//_____________________________________________________
46AliHFCorrelator::AliHFCorrelator() :
47//
48// default constructor
49//
50TNamed(),
51fPoolMgr(0x0),
52fPool(0x0),
53fhadcuts(0x0),
54fAODEvent(0x0),
55fAssociatedTracks(0x0),
56fmcArray(0x0),
57fReducedPart(0x0),
58fD0cand(0x0),
59fhypD0(0),
c84dbedf 60fDCharge(0),
bce70c96 61
62fmixing(kFALSE),
63fmontecarlo(kFALSE),
64fsystem(kFALSE),
65fselect(0),
66fUseImpactParameter(0),
67fPIDmode(0),
68
69fNofTracks(0),
70fPoolContent(0),
71
72fPhiMin(0),
73fPhiMax(0),
74
75fPtTrigger(0),
76fPhiTrigger(0),
77fEtaTrigger(0),
78
79
80fDeltaPhi(0),
81fDeltaEta(0),
82fk0InvMass(0)
83
84{
85 // default constructor
86}
87
88
89
90//_____________________________________________________
91AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :
92TNamed(name,"title"),
93fPoolMgr(0x0),
94fPool(0x0),
95fhadcuts(0x0),
96fAODEvent(0x0),
97fAssociatedTracks(0x0),
98fmcArray(0x0),
99fReducedPart(0x0),
100fD0cand(0x0),
101fhypD0(0),
c84dbedf 102fDCharge(0),
bce70c96 103
104fmixing(kFALSE),
105fmontecarlo(kFALSE),
106fsystem(ppOrPbPb),
107fselect(0),
108fUseImpactParameter(0),
109fPIDmode(0),
110
111fNofTracks(0),
112fPoolContent(0),
113
114fPhiMin(0),
115fPhiMax(0),
116
117fPtTrigger(0),
118fPhiTrigger(0),
119fEtaTrigger(0),
120
121
122fDeltaPhi(0),
123fDeltaEta(0),
124fk0InvMass(0)
125{
126 fhadcuts = cuts;
127}
128
129//_____________________________________________________
130AliHFCorrelator::~AliHFCorrelator()
131{
132//
133// destructor
134//
135
136 if(fPoolMgr) {delete fPoolMgr; fPoolMgr=0;}
137 if(fPool) {delete fPool; fPool=0;}
138 if(fhadcuts) {delete fhadcuts; fhadcuts=0;}
139 if(fAODEvent) {delete fAODEvent; fAODEvent=0;}
140 if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}
141 if(fmcArray) {delete fmcArray; fmcArray=0;}
142 if(fReducedPart) {delete fReducedPart; fReducedPart=0;}
143 if(fD0cand) {delete fD0cand; fD0cand=0;}
144
145
146 if(fNofTracks) fNofTracks = 0;
147
148 if(fPhiMin) fPhiMin = 0;
149 if(fPhiMax) fPhiMax = 0;
150
151 if(fPtTrigger) fPtTrigger=0;
152 if(fPhiTrigger) fPhiTrigger=0;
153 if(fEtaTrigger) fEtaTrigger=0;
154
155 if(fDeltaPhi) fDeltaPhi=0;
156 if(fDeltaEta) fDeltaEta=0;
157
158 if(fk0InvMass) fk0InvMass=0;
159
160}
161//_____________________________________________________
162Bool_t AliHFCorrelator::DefineEventPool(){
163 // definition of the Pool Manager for Event Mixing
164
165
166 Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();
167 Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();
168 Int_t NofCentBins = fhadcuts->GetNCentPoolBins();
169 Double_t * CentBins = fhadcuts->GetCentPoolBins();
170 Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();
171 Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();
172
173
174 fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);
175 if(!fPoolMgr) return kFALSE;
176 return kTRUE;
177}
178//_____________________________________________________
179Bool_t AliHFCorrelator::Initialize(){
180
181 cout << "AliHFCorrelator::Initialize"<< endl;
182 if(!fAODEvent) cout << "No AOD event" << endl;
183
184 AliCentrality *centralityObj = 0;
185 Int_t multiplicity = -1;
186 Double_t MultipOrCent = -1;
187
188 // initialize the pool for event mixing
189 if(!fsystem){ // pp
190 multiplicity = fAODEvent->GetNTracks();
191 MultipOrCent = multiplicity; // convert from Int_t to Double_t
192 }
193 if(fsystem){ // PbPb
194
195 centralityObj = fAODEvent->GetHeader()->GetCentralityP();
196 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
197 AliInfo(Form("Centrality is %f", MultipOrCent));
198 }
199
200 AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();
201 Double_t zvertex = vtx->GetZ(); // zvertex
202 Double_t * CentBins = fhadcuts->GetCentPoolBins();
203 Double_t poolmin=CentBins[0];
204 Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];
c84dbedf 205
bce70c96 206
207 if(TMath::Abs(zvertex)>=10 || MultipOrCent>poolmax || MultipOrCent < poolmin) {
208 if(!fsystem)AliInfo(Form("pp Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
c84dbedf 209 if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
bce70c96 210
211 return kFALSE;
212 }
213
214 fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);
215
216 if (!fPool){
217 AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));
218 return kFALSE;
219 }
c84dbedf 220 //fPool->PrintInfo();
bce70c96 221 return kTRUE;
222}
223
224//_____________________________________________________
225Bool_t AliHFCorrelator::ProcessEventPool(){
226 // analysis on Mixed Events
c84dbedf 227 //cout << "AliHFCorrelator::ProcessEventPool"<< endl;
bce70c96 228 if(!fmixing) return kFALSE;
229 if(!fPool->IsReady()) return kFALSE;
230 if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
c84dbedf 231 // fPool->PrintInfo();
bce70c96 232 fPoolContent = fPool->GetCurrentNEvents();
233
234 return kTRUE;
235
236}
237
238//_____________________________________________________
239Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex){
240
241 fAssociatedTracks = new TObjArray();
242
243 if(!fmixing){ // analysis on Single Event
244
245
246 if(fselect==1 || fselect ==2) fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
247 if(fselect==3) {fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);}
248
249
250 }
251
252 if(fmixing) { // analysis on Mixed Events
253
254
255 fAssociatedTracks = fPool->GetEvent(EventLoopIndex);
256
257
258
259
260 } // end if mixing
261
262 if(!fAssociatedTracks) return kFALSE;
263
264 fNofTracks = fAssociatedTracks->GetEntriesFast();
265
266 return kTRUE;
267
268}
269//_____________________________________________________
270Bool_t AliHFCorrelator::Correlate(Int_t loopindex){
271
272 if(loopindex >= fNofTracks) return kFALSE;
273 if(!fAssociatedTracks) return kFALSE;
274
275 fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);
276
277
278 fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());
279
280 fDeltaEta = fEtaTrigger - fReducedPart->Eta();
281
282 return kTRUE;
283
284}
285
286//_____________________________________________________
287Bool_t AliHFCorrelator::PoolUpdate(){
288
289 if(!fmixing) return kFALSE;
290 if(!fPool) return kFALSE;
291 if(fmixing) { // update the pool for Event Mixing
c84dbedf 292 TObjArray* objArr = NULL;
293 if(fselect==1 || fselect==2) objArr = (TObjArray*)AcceptAndReduceTracks(fAODEvent);
294 if(fselect==3) objArr = (TObjArray*)AcceptAndReduceKZero(fAODEvent);
295 if(objArr->GetEntriesFast()>0) fPool->UpdatePool(objArr); // updating the pool only if there are entries in the array
bce70c96 296 }
c84dbedf 297
bce70c96 298 return kTRUE;
c84dbedf 299
bce70c96 300}
301
302//_____________________________________________________
303Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){
304 Double_t pi = TMath::Pi();
305
306 if(phi<fPhiMin) phi = phi + 2*pi;
307 if(phi>fPhiMax) phi = phi - 2*pi;
308
309 return phi;
310}
311
312//_____________________________________________________
313TObjArray* AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){
314
315 Int_t nTracks = inputEvent->GetNTracks();
316 AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
317 Double_t Bz = inputEvent->GetMagneticField();
318
319
320 TObjArray* tracksClone = new TObjArray;
321 tracksClone->SetOwner(kTRUE);
322 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
323
324 AliAODTrack* track = inputEvent->GetTrack(iTrack);
325 if (!track) continue;
326 if(!fhadcuts->IsHadronSelected(track)) continue; // apply ESD level selections
c84dbedf 327 if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
328
bce70c96 329 Double_t pT = track->Pt();
330
331 //compute impact parameter
332 Double_t d0z0[2],covd0z0[3];
333 Double_t d0=-999999.;
334 if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
335 else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
336
337 if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
338 if(fUseImpactParameter==2) { // use impact parameter over resolution
bec72d8c 339 if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]);
bce70c96 340 else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
341
342 }
343
c84dbedf 344
bce70c96 345
346
c84dbedf 347 if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
bce70c96 348 Bool_t rejectsoftpi = kFALSE;
349 if(fD0cand) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0);
350
351
352 if(fselect ==2){
353 if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
354 }
355
c84dbedf 356 tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge()));
bce70c96 357 }
358
359
360 return tracksClone;
361}
362
363//_____________________________________________________
364TObjArray* AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){
365
366 Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
367 TObjArray* KZeroClone = new TObjArray;
368 AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();
369 Int_t v0label = -1;
370 Int_t pdgDgK0toPipi[2] = {211,211};
371 Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();
372 const Int_t size = inputEvent->GetNumberOfV0s();
373 Int_t prevnegID[size];
374 Int_t prevposID[size];
375 for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
376 AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);
377 if(!v0) {cout << "is not a v0 at step " << iv0 << endl; continue;}
378
379 if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
380
381 // checks to avoid double counting
382 Int_t negID = -999;
383 Int_t posID = -998;
384 //Int_t a = 0;
385 prevnegID[iv0] = -997;
386 prevposID[iv0]= -996;
387 negID = v0->GetNegID();
388 posID = v0->GetPosID();
389 Bool_t DoubleCounting = kFALSE;
390
391 for(Int_t k=0; k<iv0; k++){
392 if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){
393 DoubleCounting = kTRUE;
394 //a=k;
395 break;
396 }//end if
397 } // end for
398
c84dbedf 399 if(DoubleCounting) continue;
400 else{
401 prevposID[iv0] = posID;
402 prevnegID[iv0] = negID;
403 }
bce70c96 404
405 if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
406 Double_t k0pt = v0->Pt();
407 Double_t k0eta = v0->Eta();
408 Double_t k0Phi = v0->Phi();
409 fk0InvMass = v0->MassK0Short();
410
411 //if (loopindex == 0) {
412 // if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
413 // if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*
414 //}
415 // if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
416
417 if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
418 KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));
419
420 }
421
422 return KZeroClone;
423}