PWGJE/FlavourJetTasks/AliAnalysisTaskEmcalJetV0CF.cxx

   1 #include <THnSparse.h>
   2 #include <TClonesArray.h>
   3 #include <TParticle.h>
   4
   5 #include "AliAODEvent.h"
   6 #include "AliESDEvent.h"
   7 #include "AliMCEvent.h"
   8 #include "AliStack.h"
   9 #include "AliMCParticle.h"
  10 #include "AliAODMCParticle.h"
  11
  12 #include "AliEmcalJet.h"
  13 #include "AliJetContainer.h"
  14 #include "AliParticleContainer.h"
  15 #include "AliClusterContainer.h"
  16
  17 #include "AliPicoV0MC.h"
  18 #include "AliAnalysisTaskEmcalJetV0CF.h"
  19
  20 ClassImp(AliAnalysisTaskEmcalJetV0CF)
  21
  22 //_____________________________________________________________________________
  23 AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF() :
  24 AliAnalysisTaskEmcalJet(),
  25 fKaCutNS(0.),
  26 fLaCutNS(0.),
  27 fV0CutMinEta(0.),
  28 fV0CutMaxEta(0.),
  29 fEventAOD(0),
  30 fEventESD(0),
  31 fCentInfo(0),
  32 fJetsContRD(0),
  33 fTracksContRD(0),
  34 fCaloClustersContRD(0),
  35 fJetsContMC(0),
  36 fTracksContMC(0),
  37 fV0s(0),
  38 fHistoKshortInvM(0),
  39 fHistoLambdaInvM(0),
  40 fHistoAntiLaInvM(0),
  41 fListUserOutputs(0)
  42 {
  43 //
  44 //  AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
  45 //
  46 }
  47
  48 //_____________________________________________________________________________
  49 AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF(const char *name, Bool_t bHistos) :
  50 AliAnalysisTaskEmcalJet(name,bHistos),
  51 fKaCutNS(6.),
  52 fLaCutNS(6.),
  53 fV0CutMinEta(-10.),
  54 fV0CutMaxEta(10.),
  55 fEventAOD(0),
  56 fEventESD(0),
  57 fCentInfo(0),
  58 fJetsContRD(0),
  59 fTracksContRD(0),
  60 fCaloClustersContRD(0),
  61 fJetsContMC(0),
  62 fTracksContMC(0),
  63 fV0s(0),
  64 fHistoKshortInvM(0),
  65 fHistoLambdaInvM(0),
  66 fHistoAntiLaInvM(0),
  67 fListUserOutputs(0)
  68 {
  69 //
  70 //  AliAnalysisTaskEmcalJetV0CF::AliAnalysisTaskEmcalJetV0CF
  71 //
  72
  73   AliAnalysisTaskEmcal::fGeneralHistograms = bHistos;
  74
  75   DefineOutput(2, TList::Class());
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF()
  80 {
  81 //
  82 //  AliAnalysisTaskEmcalJetV0CF::~AliAnalysisTaskEmcalJetV0CF
  83 //
  84
  85   if (fEventAOD) { delete fEventAOD; fEventAOD = 0; }
  86   if (fEventESD) { delete fEventESD; fEventESD = 0; }
  87   if (fCentInfo) { delete fCentInfo; fCentInfo = 0; }
  88
  89   if (fJetsContRD)         { delete fJetsContRD;         fJetsContRD         = 0; }
  90   if (fTracksContRD)       { delete fTracksContRD;       fTracksContRD       = 0; }
  91   if (fCaloClustersContRD) { delete fCaloClustersContRD; fCaloClustersContRD = 0; }
  92
  93   if (fJetsContMC)   { delete fJetsContMC;   fJetsContMC   = 0; }
  94   if (fTracksContMC) { delete fTracksContMC; fTracksContMC = 0; }
  95
  96   if (fV0s) { delete fV0s; fV0s = 0; }
  97
  98   if (fHistoKshortInvM) { delete fHistoKshortInvM; fHistoKshortInvM = 0; }
  99   if (fHistoLambdaInvM) { delete fHistoLambdaInvM; fHistoLambdaInvM = 0; }
 100   if (fHistoAntiLaInvM) { delete fHistoAntiLaInvM; fHistoAntiLaInvM = 0; }
 101
 102   if (fListUserOutputs) { delete fListUserOutputs; fListUserOutputs = 0; }
 103 }
 104
 105 //_____________________________________________________________________________
 106 void AliAnalysisTaskEmcalJetV0CF::Init()
 107 {
 108 //
 109 //  AliAnalysisTaskEmcalJetV0CF::Init
 110 //
 111
 112   return;
 113 }
 114
 115 //_____________________________________________________________________________
 116 void AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects()
 117 {
 118 //
 119 //  AliAnalysisTaskEmcalJetV0CF::UserCreateOutputObjects
 120 //
 121
 122   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 123 //=============================================================================
 124
 125   fJetsContRD = GetJetContainer(0);
 126   if (fJetsContRD) {
 127     fTracksContRD       = fJetsContRD->GetParticleContainer();
 128     fCaloClustersContRD = fJetsContRD->GetClusterContainer();
 129   }
 130
 131   fJetsContMC = GetJetContainer(1);
 132   if (fJetsContMC) fTracksContMC = fJetsContMC->GetParticleContainer();
 133 //=============================================================================
 134
 135   fListUserOutputs = new TList();
 136   fListUserOutputs->SetOwner();
 137   CreateUserOutputHistograms();
 138   PostData(2, fListUserOutputs);
 139   return;
 140 }
 141
 142 //_____________________________________________________________________________
 143 void AliAnalysisTaskEmcalJetV0CF::Terminate(Option_t *opt)
 144 {
 145 //
 146 //  AliAnalysisTaskEmcalJetV0CF::Terminate
 147 //
 148
 149   AliAnalysisTaskEmcalJet::Terminate(opt);
 150
 151   return;
 152 }
 153
 154 //_____________________________________________________________________________
 155 Bool_t AliAnalysisTaskEmcalJetV0CF::Run()
 156 {
 157 //
 158 //  AliAnalysisTaskEmcalJetV0CF::Run
 159 //
 160
 161   if (!AliAnalysisTaskEmcalJet::Run()) return kFALSE;
 162
 163   return kTRUE;
 164 }
 165
 166 //_____________________________________________________________________________
 167 Bool_t AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects()
 168 {
 169 //
 170 //  AliAnalysisTaskEmcalJetV0CF::RetrieveEventObjects
 171 //
 172
 173   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
 174
 175   fCentInfo = InputEvent()->GetCentrality(); if (!fCentInfo) return kFALSE;
 176
 177   fEventAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 178   fEventESD = dynamic_cast<AliESDEvent*>(InputEvent());
 179   if ((!fEventAOD) && (!fEventESD)) return kFALSE;
 180
 181   return kTRUE;
 182 }
 183
 184 //_____________________________________________________________________________
 185 Bool_t AliAnalysisTaskEmcalJetV0CF::IsEventSelected()
 186 {
 187 //
 188 //  AliAnalysisTaskEmcalJetV0CF::IsEventSelected
 189 //
 190
 191   if (!AliAnalysisTaskEmcalJet::IsEventSelected()) return kFALSE;
 192
 193   return kTRUE;
 194 }
 195
 196 //_____________________________________________________________________________
 197 Bool_t AliAnalysisTaskEmcalJetV0CF::FillHistograms()
 198 {
 199 //
 200 //  AliAnalysisTaskEmcalJetV0CF::FillHistograms
 201 //
 202
 203   if (!AliAnalysisTaskEmcalJet::FillHistograms()) return kFALSE;
 204
 205   if (FillRecoInfo()) return kFALSE;
 206   if (FillKineInfo()) return kFALSE;
 207
 208   return kTRUE;
 209 }
 210
 211 //_____________________________________________________________________________
 212 Bool_t AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms()
 213 {
 214 //
 215 //  AliAnalysisTaskEmcalJetV0CF::FillGeneralHistograms
 216 //
 217
 218   if (!AliAnalysisTaskEmcalJet::FillGeneralHistograms()) return kFALSE;
 219
 220   return kTRUE;
 221 }
 222
 223 //_____________________________________________________________________________
 224 void AliAnalysisTaskEmcalJetV0CF::ExecOnce()
 225 {
 226 //
 227 //  AliAnalysisTaskEmcalJetV0CF::ExecOnce
 228 //
 229
 230   AliAnalysisTaskEmcalJet::ExecOnce();
 231
 232   if (!fInitialized) return;
 233
 234   if (!fV0s) {
 235     fV0s = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoV0s"));
 236
 237     if (!fV0s) {
 238       AliError(Form("%s: Could not retrieve V0 %s!", GetName(), "PicoV0s"));
 239       fInitialized = kFALSE;
 240       return;
 241     }
 242   }
 243
 244   return;
 245 }
 246
 247 //_____________________________________________________________________________
 248 void AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms()
 249 {
 250 //
 251 //  AliAnalysisTaskEmcalJetV0CF::CreateUserOutputHistograms
 252 //
 253
 254   if (!fListUserOutputs) return;
 255
 256   Bool_t bStatusTmpH = TH1::AddDirectoryStatus();
 257   TH1::AddDirectory(kFALSE);
 258
 259   const Int_t nV0 = 8; // 0: particle type
 260                        //    ==0, Kshort
 261                        //    ==1, Lambda
 262                        //    ==2, AntiLa
 263                        // 1: Jet pT bin
 264                        //    == 0.5, jet pT>10.
 265                        //    == 1.5, jet pT>15.
 266                        //    == 2.5, jet pT>20.
 267                        //    == 3.5, jet pT>25.
 268                        // 2: V0M
 269                        // 3: V0A
 270                        // 4: CL1
 271                        // 5: ZNA
 272                        // 6: eta
 273                        // 7: Pt
 274   const Int_t    nV0Bin[nV0] = {  3,   4,  210,  210,  210,  210, 100, 1000  };
 275   const Double_t dV0Min[nV0] = { -0.5, 0., -10., -10., -10., -10., -5.,   0. };
 276   const Double_t dV0Max[nV0] = {  2.5, 4., 200., 200., 200., 200.,  5., 100. };
 277
 278   THnSparseD *hsReco = new THnSparseD("hsReco", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsReco);
 279   THnSparseD *hsKine = new THnSparseD("hsKine", "", nV0, nV0Bin, dV0Min, dV0Max); fListUserOutputs->Add(hsKine);
 280
 281   TH1::AddDirectory(bStatusTmpH);
 282   return;
 283 }
 284
 285 //_____________________________________________________________________________
 286 Bool_t AliAnalysisTaskEmcalJetV0CF::FillRecoInfo()
 287 {
 288 //
 289 //  AliAnalysisTaskEmcalJetV0CF::FillRecoInfo
 290 //
 291
 292   if (!fV0s) return kTRUE;
 293 //=============================================================================
 294
 295   Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
 296   Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
 297   Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
 298   Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");
 299
 300   THnSparseD *hs = dynamic_cast<THnSparseD*>(fListUserOutputs->FindObject("hsReco"));
 301   if (!hs)
 302     return kTRUE;  // should not happen; make Coverity happen
 303 //=============================================================================
 304
 305   AliPicoV0MC *pV0 = 0;
 306   for (Int_t i=0; i<fV0s->GetEntriesFast(); i++) {
 307     pV0 = static_cast<AliPicoV0MC*>(fV0s->At(i)); if (!pV0) continue;
 308
 309     if (!pV0->IsV0PhysicalPrimary()) { pV0 = 0; continue; }
 310     if (!pV0->IsV0InEtaAcc(fV0CutMinEta,fV0CutMaxEta)) { pV0 = 0; continue; }
 311
 312     TVector3 vV0 = pV0->KineMC().Vect();
 313     Double_t dPt = vV0.Pt();
 314     Double_t dVar[8];
 315
 316     dVar[0] = -1.;
 317     if (pV0->IsKshort()) {
 318       Int_t k = fHistoKshortInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
 319
 320       Double_t dMean  = fHistoKshortInvM->GetBinContent(k);
 321       Double_t dSigma = fHistoKshortInvM->GetBinError(k);
 322
 323       Double_t dUpperL = dMean - (fKaCutNS * dSigma);
 324       Double_t dLowerR = dMean + (fKaCutNS * dSigma);
 325
 326       Double_t dInvM = pV0->KineKshort().M();
 327       if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }
 328
 329       dVar[0] = 0.;
 330     }
 331
 332     if (pV0->IsLambda()) {
 333       Int_t k = fHistoLambdaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
 334
 335       Double_t dMean  = fHistoLambdaInvM->GetBinContent(k);
 336       Double_t dSigma = fHistoLambdaInvM->GetBinError(k);
 337
 338       Double_t dUpperL = dMean - (fLaCutNS * dSigma);
 339       Double_t dLowerR = dMean + (fLaCutNS * dSigma);
 340
 341       Double_t dInvM = pV0->KineLambda().M();
 342       if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }
 343
 344       dVar[0] = 1.;
 345     }
 346
 347     if (pV0->IsAntiLa()) {
 348       Int_t k = fHistoAntiLaInvM->FindBin(dPt); if (k<=0) { pV0 = 0; continue; }
 349
 350       Double_t dMean  = fHistoAntiLaInvM->GetBinContent(k);
 351       Double_t dSigma = fHistoAntiLaInvM->GetBinError(k);
 352
 353       Double_t dUpperL = dMean - (fLaCutNS * dSigma);
 354       Double_t dLowerR = dMean + (fLaCutNS * dSigma);
 355
 356       Double_t dInvM = pV0->KineAntiLa().M();
 357       if ((dInvM<dUpperL) || (dInvM>=dLowerR)) { pV0 = 0; continue; }
 358
 359       dVar[0] = 2.;
 360     }
 361
 362     if (dVar[0]<-0.5) { pV0 = 0; continue; }
 363
 364     dVar[1] = -1.;
 365     if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
 366     if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
 367     if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
 368     if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
 369     if (dVar[1]<0.) { pV0 = 0; continue; }
 370
 371     dVar[2] = dV0M;
 372     dVar[3] = dV0A;
 373     dVar[4] = dCL1;
 374     dVar[5] = dZNA;
 375
 376     dVar[6] = vV0.Eta();
 377     dVar[7] = dPt;
 378
 379     hs->Fill(dVar);
 380
 381     pV0 = 0;
 382   }
 383
 384   return kFALSE;
 385 }
 386
 387 //_____________________________________________________________________________
 388 Bool_t AliAnalysisTaskEmcalJetV0CF::FillKineInfo()
 389 {
 390 //
 391 //  AliAnalysisTaskEmcalJetV0CF::FillKineInfo
 392 //
 393
 394   Double_t dV0M = fCentInfo->GetCentralityPercentile("V0M");
 395   Double_t dV0A = fCentInfo->GetCentralityPercentile("V0A");
 396   Double_t dCL1 = fCentInfo->GetCentralityPercentile("CL1");
 397   Double_t dZNA = fCentInfo->GetCentralityPercentile("ZNA");
 398
 399   THnSparseD *hs = dynamic_cast<THnSparseD*>(fListUserOutputs->FindObject("hsKine"));
 400   if (!hs) {
 401     return kTRUE; // Should not happen; make Coverity happy
 402   }
 403 //=============================================================================
 404
 405   AliStack *pStack = 0;
 406   if (fEventESD) { pStack = MCEvent()->Stack(); if (!pStack) return kTRUE; }
 407 //=============================================================================
 408
 409   TParticle        *pESD = 0;
 410   AliAODMCParticle *pAOD = 0;
 411   for (Int_t i=0; i<MCEvent()->GetNumberOfTracks(); i++) {
 412     if (fEventAOD) { pAOD = (AliAODMCParticle*)MCEvent()->GetTrack(i);              if (!pAOD) continue; }
 413     if (fEventESD) { pESD =   ((AliMCParticle*)MCEvent()->GetTrack(i))->Particle(); if (!pESD) continue; }
 414
 415     Bool_t bPhy = kFALSE;
 416     if (pAOD) bPhy =   pAOD->IsPhysicalPrimary();
 417     if (pESD) bPhy = pStack->IsPhysicalPrimary(i);
 418     if (!bPhy) { pAOD=0; pESD=0; continue; }
 419
 420     Double_t  dEta = 0.;
 421     if (pAOD) dEta = pAOD->Eta();
 422     if (pESD) dEta = pESD->Eta();
 423     if ((dEta<fV0CutMinEta) || (dEta>=fV0CutMaxEta)) { pAOD=0; pESD=0; continue; }
 424
 425     Int_t id = 0;
 426     if (pAOD) id = pAOD->GetPdgCode();
 427     if (pESD) id = pESD->GetPdgCode();
 428
 429     Double_t dVar[8]; dVar[0] = -1.;
 430     if (id== 310 ) dVar[0] = 0.;
 431     if (id== 3122) dVar[0] = 1.;
 432     if (id==-3122) dVar[0] = 2.;
 433     if (dVar[0]<-0.5) { pAOD=0; pESD=0; continue; }
 434
 435     TVector3 vV0;
 436     if (pAOD) vV0.SetXYZ(pAOD->Px(), pAOD->Py(), pAOD->Pz());
 437     if (pESD) vV0.SetXYZ(pESD->Px(), pESD->Py(), pESD->Pz());
 438
 439     dVar[1] = -1.;
 440     if (IsV0InJet(vV0,10.)) dVar[1] = 0.5;
 441     if (IsV0InJet(vV0,15.)) dVar[1] = 1.5;
 442     if (IsV0InJet(vV0,20.)) dVar[1] = 2.5;
 443     if (IsV0InJet(vV0,25.)) dVar[1] = 3.5;
 444     if (dVar[1]<0.) { pAOD=0; pESD=0; continue; }
 445
 446     dVar[2] = dV0M;
 447     dVar[3] = dV0A;
 448     dVar[4] = dCL1;
 449     dVar[5] = dZNA;
 450     dVar[6] = dEta;
 451     if (pAOD) dVar[7] = pAOD->Pt();
 452     if (pESD) dVar[7] = pESD->Pt();
 453     hs->Fill(dVar);
 454
 455     pAOD = 0;
 456     pESD = 0;
 457   }
 458
 459   return kFALSE;
 460 }
 461
 462 //_____________________________________________________________________________
 463 Bool_t AliAnalysisTaskEmcalJetV0CF::IsV0InJet(TVector3 vV0, Double_t dJetPtMin)
 464 {
 465 //
 466 //  AliAnalysisTaskEmcalJetV0CF::IsV0InJet
 467 //
 468
 469   if (!fJetsContRD) return kFALSE;
 470 //=============================================================================
 471
 472   TVector3 vJet;
 473   Double_t dJetRadius = fJetsContRD->GetJetRadius();
 474   AliEmcalJet *pJet = fJetsContRD->GetNextAcceptJet(0); while (pJet) {
 475     Double_t dPt = fJetsContRD->GetJetPtCorr(fJetsContRD->GetCurrentID());
 476     if (dPt<dJetPtMin) { pJet = fJetsContRD->GetNextAcceptJet(); continue; }
 477
 478     vJet.SetPtEtaPhi(dPt, pJet->Eta(), pJet->Phi());
 479     if (vJet.DeltaR(vV0)<dJetRadius) return kTRUE;
 480     pJet = fJetsContRD->GetNextAcceptJet();
 481   }
 482
 483   return kFALSE;
 484 }