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