Fixes for std:: need with the trunk of Root (Jochen, Yves)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliHFCorrelator.cxx
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
42 using std::cout;
43 using std::endl;
44
45 //_____________________________________________________
46 AliHFCorrelator::AliHFCorrelator() :
47 //
48 // default constructor
49 //
50 TNamed(),
51 fPoolMgr(0x0),         
52 fPool(0x0),
53 fhadcuts(0x0),
54 fAODEvent(0x0),
55 fAssociatedTracks(0x0),
56 fmcArray(0x0),
57 fReducedPart(0x0),
58 fD0cand(0x0), 
59 fhypD0(0), 
60 fDCharge(0),
61
62 fmixing(kFALSE),
63 fmontecarlo(kFALSE),
64 fsystem(kFALSE),
65 fselect(0),
66 fUseImpactParameter(0),
67 fPIDmode(0),
68
69 fNofTracks(0),
70 fPoolContent(0),
71
72 fPhiMin(0),
73 fPhiMax(0),
74
75 fPtTrigger(0),
76 fPhiTrigger(0),
77 fEtaTrigger(0),
78
79
80 fDeltaPhi(0),
81 fDeltaEta(0),
82 fk0InvMass(0)
83
84 {
85         // default constructor  
86 }
87
88
89
90 //_____________________________________________________
91 AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :
92 TNamed(name,"title"),
93 fPoolMgr(0x0),         
94 fPool(0x0),
95 fhadcuts(0x0),
96 fAODEvent(0x0),
97 fAssociatedTracks(0x0),
98 fmcArray(0x0),
99 fReducedPart(0x0),
100 fD0cand(0x0), 
101 fhypD0(0),
102 fDCharge(0),
103
104 fmixing(kFALSE),
105 fmontecarlo(kFALSE),
106 fsystem(ppOrPbPb),
107 fselect(0),
108 fUseImpactParameter(0),
109 fPIDmode(0),
110
111 fNofTracks(0),
112 fPoolContent(0),
113
114 fPhiMin(0),
115 fPhiMax(0),
116
117 fPtTrigger(0),
118 fPhiTrigger(0),
119 fEtaTrigger(0),
120
121
122 fDeltaPhi(0),
123 fDeltaEta(0),
124 fk0InvMass(0)
125 {
126         fhadcuts = cuts; 
127 }
128
129 //_____________________________________________________
130 AliHFCorrelator::~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 //_____________________________________________________
162 Bool_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 //_____________________________________________________
179 Bool_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()];
205
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));
209                 if(fsystem) AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f  out of pool bounds, SKIPPING",zvertex,MultipOrCent));
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         }
220         //fPool->PrintInfo();
221         return kTRUE;
222 }
223
224 //_____________________________________________________
225 Bool_t AliHFCorrelator::ProcessEventPool(){
226          // analysis on Mixed Events
227         //cout << "AliHFCorrelator::ProcessEventPool"<< endl;
228                 if(!fmixing) return kFALSE;
229                 if(!fPool->IsReady()) return kFALSE;
230                 if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
231         //      fPool->PrintInfo();
232                 fPoolContent = fPool->GetCurrentNEvents();
233                 
234                 return kTRUE;
235         
236 }
237
238 //_____________________________________________________
239 Bool_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 //_____________________________________________________
270 Bool_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 //_____________________________________________________
287 Bool_t AliHFCorrelator::PoolUpdate(){
288
289         if(!fmixing) return kFALSE;
290         if(!fPool) return kFALSE;
291         if(fmixing) { // update the pool for Event Mixing
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
296         }
297                 
298         return kTRUE;
299         
300 }
301                 
302 //_____________________________________________________
303 Double_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 //_____________________________________________________
313 TObjArray*  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
327         if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
328                 
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
339                   if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]); 
340                   else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
341                         
342                 }
343                 
344         
345                 
346                 
347                 if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
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                 
356                 tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge()));
357         }
358         
359         
360         return tracksClone;
361 }
362
363 //_____________________________________________________
364 TObjArray*  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                 
399                 if(DoubleCounting) continue;
400                 else{ 
401                         prevposID[iv0] = posID; 
402                         prevnegID[iv0] = negID;
403                 }
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 }