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