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