PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.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 DStar - Hadron Correlations Analysis
  18 //
  19 //-----------------------------------------------------------------------
  20 //
  21 //
  22 //                                                 Author S.Bjelogrlic
  23 //                         Utrecht University
  24 //                      sandro.bjelogrlic@cern.ch
  25 //
  26 //-----------------------------------------------------------------------
  27
  28 /* $Id: AliAnalysisTaskDStarCorrelations.cxx 65139 2013-11-25 14:47:45Z fprino $ */
  29
  30 //#include <TDatabasePDG.h>
  31 #include <TParticle.h>
  32 #include <TVector3.h>
  33 #include <TChain.h>
  34 #include "TROOT.h"
  35
  36 #include "AliAnalysisTaskDStarCorrelations.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliHFAssociatedTrackCuts.h"
  39 #include "AliAODRecoDecay.h"
  40 #include "AliAODRecoCascadeHF.h"
  41 #include "AliAODRecoDecayHF2Prong.h"
  42 #include "AliAODPidHF.h"
  43 #include "AliVParticle.h"
  44 #include "AliAnalysisManager.h"
  45 #include "AliAODInputHandler.h"
  46 #include "AliAODHandler.h"
  47 #include "AliESDtrack.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliNormalizationCounter.h"
  50 #include "AliReducedParticle.h"
  51 #include "AliHFCorrelator.h"
  52 #include "AliAODMCHeader.h"
  53 #include "AliEventPoolManager.h"
  54 #include "AliVertexingHFUtils.h"
  55
  56 using std::cout;
  57 using std::endl;
  58
  59
  60 ClassImp(AliAnalysisTaskDStarCorrelations)
  61
  62
  63 //__________________________________________________________________________
  64 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
  65 AliAnalysisTaskSE(),
  66 fhandler(0x0),
  67 fmcArray(0x0),
  68 fCounter(0x0),
  69 fCorrelator(0x0),
  70 fselect(0),
  71 fmontecarlo(kFALSE),
  72 fmixing(kFALSE),
  73 fFullmode(kFALSE),
  74 fSystem(pp),
  75 fEfficiencyVariable(kNone),
  76 fBkgMethod(kDZeroSB),
  77 fReco(kTRUE),
  78 fUseEfficiencyCorrection(kFALSE),
  79 fUseDmesonEfficiencyCorrection(kFALSE),
  80 fUseCentrality(kFALSE),
  81 fUseHadronicChannelAtKineLevel(kFALSE),
  82 fPhiBins(32),
  83 fEvents(0),
  84 fDebugLevel(0),
  85 fDisplacement(0),
  86 fDim(0),
  87 fNofPtBins(0),
  88 fMaxEtaDStar(0.9),
  89 fDMesonSigmas(0),
  90
  91 fD0Window(0),
  92
  93 fOutput(0x0),
  94 fDmesonOutput(0x0),
  95 fTracksOutput(0x0),
  96 fEMOutput(0x0),
  97 fCorrelationOutput(0x0),
  98 fOutputMC(0x0),
  99
 100 fCuts(0),
 101 fAssocCuts(0),
 102 fUtils(0),
 103 fTracklets(0),
 104 fDeffMapvsPt(0),
 105 fDeffMapvsPtvsMult(0),
 106 fDeffMapvsPtvsEta(0)
 107 {
 108     SetDim();
 109     // default constructor
 110 }
 111
 112 //__________________________________________________________________________
 113 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :
 114 AliAnalysisTaskSE(name),
 115
 116 fhandler(0x0),
 117 fmcArray(0x0),
 118 fCounter(0x0),
 119 fCorrelator(0x0),
 120 fselect(0),
 121 fmontecarlo(kFALSE),
 122 fmixing(kFALSE),
 123 fFullmode(mode),
 124 fSystem(syst),
 125 fEfficiencyVariable(kNone),
 126 fBkgMethod(kDZeroSB),
 127 fReco(kTRUE),
 128 fUseEfficiencyCorrection(kFALSE),
 129 fUseDmesonEfficiencyCorrection(kFALSE),
 130 fUseCentrality(kFALSE),
 131 fUseHadronicChannelAtKineLevel(kFALSE),
 132 fPhiBins(32),
 133 fEvents(0),
 134 fDebugLevel(0),
 135 fDisplacement(0),
 136 fDim(0),
 137 fNofPtBins(0),
 138 fMaxEtaDStar(0.9),
 139 fDMesonSigmas(0),
 140 fD0Window(0),
 141
 142 fOutput(0x0),
 143 fDmesonOutput(0x0),
 144 fTracksOutput(0x0),
 145 fEMOutput(0x0),
 146 fCorrelationOutput(0x0),
 147 fOutputMC(0x0),
 148
 149 fCuts(0),
 150 fAssocCuts(AsscCuts),
 151 fUtils(0),
 152 fTracklets(0),
 153 fDeffMapvsPt(0),
 154 fDeffMapvsPtvsMult(0),
 155 fDeffMapvsPtvsEta(0)
 156 {
 157      Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
 158
 159     SetDim();
 160     if(fSystem == pp)  fUseCentrality = kFALSE; else fUseCentrality = kTRUE;
 161
 162     if(fSystem == AA) fBkgMethod = kDStarSB; else fBkgMethod = kDZeroSB;
 163
 164     fCuts=cuts;
 165     fNofPtBins= fCuts->GetNPtBins();
 166     //cout << "Enlarging the DZero window " << endl;
 167     EnlargeDZeroMassWindow();
 168    // cout << "Done" << endl;
 169
 170
 171   DefineInput(0, TChain::Class());
 172
 173   DefineOutput(1,TList::Class()); // histos from data and MC
 174   DefineOutput(2,TList::Class()); // histos from MC
 175     DefineOutput(3,TList::Class()); // histos from data and MC
 176     DefineOutput(4,TList::Class()); // histos from MC
 177     DefineOutput(5,TList::Class()); // histos from data and MC
 178       DefineOutput(6,TList::Class()); // histos from data and MC
 179     DefineOutput(7,AliNormalizationCounter::Class());   // normalization
 180
 181   DefineOutput(8,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
 182   DefineOutput(9,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
 183 }
 184
 185 //__________________________________________________________________________
 186
 187 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
 188   //
 189         // destructor
 190         //
 191
 192         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");
 193
 194         if(fhandler) {delete fhandler; fhandler = 0;}
 195         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}
 196         if(fmcArray) {delete fmcArray; fmcArray = 0;}
 197         if(fCounter) {delete fCounter; fCounter = 0;}
 198         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
 199         if(fOutput) {delete fOutput; fOutput = 0;}
 200         if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
 201         if(fCuts) {delete fCuts; fCuts = 0;}
 202         if(fAssocCuts) {delete fAssocCuts; fAssocCuts=0;}
 203     if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}
 204     if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}
 205     if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}
 206
 207 }
 208
 209 //___________________________________________________________
 210 void AliAnalysisTaskDStarCorrelations::Init(){
 211         //
 212         // Initialization
 213         //
 214         if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
 215
 216         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
 217
 218
 219
 220         // Post the D* cuts
 221         PostData(8,copyfCuts);
 222
 223         // Post the hadron cuts
 224         PostData(9,fAssocCuts);
 225
 226         return;
 227 }
 228
 229
 230 //_________________________________________________
 231 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
 232         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 233
 234         //slot #1
 235         //OpenFile(0);
 236         fOutput = new TList();
 237         fOutput->SetOwner();
 238
 239         fOutputMC = new TList();
 240         fOutputMC->SetOwner();
 241
 242     fDmesonOutput = new TList();
 243         fDmesonOutput->SetOwner();
 244
 245     fTracksOutput = new TList();
 246         fTracksOutput->SetOwner();
 247
 248     fEMOutput = new TList();
 249         fEMOutput->SetOwner();
 250
 251     fCorrelationOutput = new TList();
 252         fCorrelationOutput->SetOwner();
 253
 254         // define histograms
 255         DefineHistoForAnalysis();
 256     DefineThNSparseForAnalysis();
 257
 258
 259
 260
 261         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(6)->GetContainer()->GetName()));
 262         fCounter->Init();
 263
 264     Double_t Pi = TMath::Pi();
 265     //AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t useCentrality, AliRDHFCuts *cutObject)
 266         fCorrelator = new AliHFCorrelator("Correlator",fAssocCuts,fUseCentrality,fCuts); // fAssocCuts is the hadron cut object, fSystem to switch between pp or PbPb
 267         fCorrelator->SetDeltaPhiInterval(  -0.5*Pi, 1.5*Pi); // set correct phi interval
 268         //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval
 269         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
 270         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
 271         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
 272         fCorrelator->SetUseMC(fmontecarlo);
 273         fCorrelator->SetUseReco(fReco);
 274     //  fCorrelator->SetKinkRemoval(kTRUE);
 275         Bool_t pooldef = fCorrelator->DefineEventPool();
 276
 277         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
 278
 279     fUtils = new AliAnalysisUtils();
 280
 281
 282
 283         PostData(1,fOutput); // set the outputs
 284         PostData(2,fDmesonOutput); // set the outputs
 285     PostData(3,fTracksOutput); // set the outputs
 286     PostData(4,fEMOutput); // set the outputs
 287     PostData(5,fCorrelationOutput); // set the outputs
 288     PostData(6,fOutputMC); // set the outputs
 289         PostData(7,fCounter); // set the outputs
 290 }
 291
 292 //________________________________________  user exec ____________
 293 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
 294     // cout << "Task debug check 1 " << endl;
 295      // ********************************************** LOAD THE EVENT ****************************************************
 296
 297     if (!fInputEvent) {
 298         Error("UserExec","NO EVENT FOUND!");
 299                 return;
 300     }
 301
 302     AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 303     if(!aodEvent){
 304         AliError("AOD event not found!");
 305         return;
 306     }
 307
 308     fEvents++; // event counter
 309     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
 310
 311     fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); // store event in AliNormalizationCounter
 312
 313     // load MC array
 314     fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 315     if(fmontecarlo && !fmcArray){
 316         AliError("Array of MC particles not found");
 317         return;
 318     }
 319
 320     // ********************************************** START EVENT SELECTION ****************************************************
 321
 322     Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
 323
 324     if(!isEvSel) {
 325
 326         if(fCuts->IsEventRejectedDueToPileup()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);                        cout << "Reject PILEUP" << endl;}
 327         if(fCuts->IsEventRejectedDueToCentrality()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);                    cout << "Reject CENTRALITY" << endl;}
 328         if(fCuts->IsEventRejectedDueToNotRecoVertex()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);                 cout << "Reject NOT RECO VTX" << endl;}
 329         if(fCuts->IsEventRejectedDueToVertexContributors()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);            cout << "Reject VTX CONTRIB" << endl;}
 330         if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);  cout << "Reject VTX no fid reg " << endl;}
 331         if(fCuts->IsEventRejectedDueToTrigger()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);                       cout << "Reject TRIGGER" << endl;}
 332         if(fCuts->IsEventRejectedDuePhysicsSelection()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);                cout << "Reject PHYS SEL" << endl;}
 333
 334         return;
 335     }
 336     // added event selection for pA
 337
 338     if(fSystem == pA){
 339
 340         if(fUtils->IsFirstEventInChunk(aodEvent)) {
 341             AliInfo("Rejecting the event - first in the chunk");
 342             ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);
 343             return;
 344         }
 345         if(!fUtils->IsVertexSelected2013pA(aodEvent)) {
 346             ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);
 347             AliInfo("Rejecting the event - bad vertex");
 348             return;
 349         }
 350     }
 351
 352     fCorrelator->SetAODEvent(aodEvent); // set the event in the correlator class
 353
 354     Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event within the pool definition
 355         if(!correlatorON) {
 356         AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
 357         return; // if not the case, skips the event
 358         }
 359
 360     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); // count events that are passing the event selection
 361
 362
 363
 364     // cout << "Task debug check 2 " << endl;
 365      // ********************************************** CENTRALITY/MULTIPLICITY  ****************************************************
 366
 367
 368     Double_t MultipOrCent = fCorrelator->GetCentrality(); // returns centrality or multiplicity
 369
 370
 371
 372     // ********************************************** MC SETUP  ****************************************************
 373
 374     if(fmontecarlo) {
 375
 376         fCorrelator->SetMCArray(fmcArray); // export MC array into correlatior class
 377
 378         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 379         if (fmontecarlo && !mcHeader) {
 380             AliError("Could not find MC Header in AOD");
 381             return;
 382         }
 383
 384       // select MC event type (PP, GS etc) - those are defined in the associated tracks cut object
 385         Bool_t isMCeventgood = kFALSE;
 386
 387
 388         Int_t eventType = mcHeader->GetEventType();
 389         Int_t NMCevents = fAssocCuts->GetNofMCEventType();
 390
 391         for(Int_t k=0; k<NMCevents; k++){
 392             Int_t * MCEventType = fAssocCuts->GetMCEventType();
 393
 394             if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
 395             ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
 396         }
 397
 398         if(NMCevents && !isMCeventgood){
 399             if(fDebugLevel)     std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
 400             return;
 401         }
 402
 403
 404     }// end if fmontecarlo
 405
 406
 407      // ********************************************** EVENT PROPERTIES  ****************************************************
 408     // checks on vertex and multiplicity of the event
 409         AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
 410         Double_t zVtxPosition = vtx->GetZ(); // zvertex
 411
 412
 413
 414     if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 415
 416     if(!fmixing) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);
 417     if(!fmixing) ((TH1D*)fOutput->FindObject("MultiplicityOnlyCheck"))->Fill(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
 418
 419
 420    // Int_t poolbin = fAssocCuts->GetPoolBin(MultipOrCent, zVtxPosition); // get the event pool bin - commented for the moment - to be checked
 421
 422
 423
 424      // ********************************************** D * selection  ****************************************************
 425
 426
 427     TClonesArray *arrayDStartoD0pi=0;
 428         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 429         // In case there is an AOD handler writing a standard AOD, use the AOD
 430         // event in memory rather than the input (ESD) event.
 431         aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 432         // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 433         // have to taken from the AOD event hold by the AliAODExtension
 434         AliAODHandler* aodHandler = (AliAODHandler*)
 435             ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 436         if(aodHandler->GetExtensions()) {
 437             AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 438             AliAODEvent *aodFromExt = ext->GetAOD();
 439             arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 440         }
 441         } else {
 442         arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 443         }
 444
 445
 446
 447     // D* related variables
 448
 449     Double_t ptDStar;   // D* pt
 450         Double_t phiDStar;  // D* phi
 451         Double_t etaDStar;  // D* eta
 452         Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; // flags for selection
 453         Double_t invMassDZero; // D0 inv mass
 454         Double_t deltainvMDStar; // D* - D0 invarian mass
 455
 456     Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
 457         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
 458
 459
 460         //MC tagging for DStar
 461         //D* and D0 prongs needed to MatchToMC method
 462         Int_t pdgDgDStartoD0pi[2]={421,211};
 463         Int_t pdgDgD0toKpi[2]={321,211};
 464
 465         //Bool_t isDStarCand = kFALSE;
 466     Bool_t isDfromB = kFALSE;
 467         Bool_t isEventMixingFilledPeak = kFALSE;
 468         Bool_t isEventMixingFilledSB = kFALSE;
 469     Bool_t EventHasDStarCandidate = kFALSE;
 470     Bool_t EventHasDZeroSideBandCandidate = kFALSE;
 471     Bool_t EventHasDStarSideBandCandidate = kFALSE;
 472     Bool_t isTrackForPeakFilled = kFALSE;
 473     Bool_t isTrackForSBFilled = kFALSE;
 474         //loop on D* candidates
 475
 476         Int_t looponDCands = 0;
 477         if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); // load N of D* candidates from AOD array
 478         if(!fReco) looponDCands = fmcArray->GetEntriesFast(); // load N of D* candidates from MC array
 479
 480         Int_t nOfDStarCandidates = 0;  // D candidates counter
 481         Int_t nOfSBCandidates = 0;     // SB candidates counter
 482
 483         Double_t DmesonEfficiency = 1.; // Efficiency of D meson
 484         Double_t DmesonWeight = 1.;    // Applyed weight
 485     Double_t efficiencyvariable = -999; // Variable to be used (defined by the AddTask)
 486     Double_t dmDStarWindow = 0.0019/3; // DStar window
 487
 488     // cout << "Task debug check 3 " << endl;
 489     // loop over D meson candidates
 490     for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
 491
 492         // initialize variables
 493         isInPeak = kFALSE;
 494         isInDStarSideBand = kFALSE;
 495         isInDZeroSideBand = kFALSE;
 496
 497         EventHasDStarCandidate = kFALSE;
 498         EventHasDZeroSideBandCandidate = kFALSE;
 499         EventHasDStarSideBandCandidate = kFALSE;
 500
 501
 502         isDStarMCtag = kFALSE;
 503         isDfromB = kFALSE;
 504         ptDStar = -123.4;
 505         phiDStar = -999;
 506         etaDStar = -56.;
 507         invMassDZero = - 999;
 508         deltainvMDStar = -998;
 509         AliAODRecoCascadeHF* dstarD0pi;
 510         AliAODRecoDecayHF2Prong* theD0particle;
 511         AliAODMCParticle* DStarMC=0x0;
 512         Short_t daughtercharge = -2;
 513         Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
 514         Int_t trackiddaugh1 = -1;
 515         Int_t trackidsoftPi = -1;
 516         Int_t ptbin = -1;
 517
 518         // ============================================ using reconstruction on Data or MC
 519         if(fReco){
 520             // cout << "Task debug check 4 " << endl;
 521             // get the D objects
 522             dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 523             if(!dstarD0pi->GetSecondaryVtx()) continue;
 524             theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 525             if (!theD0particle) continue;
 526
 527
 528             // apply topological cuts
 529
 530             // track quality cuts
 531             Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 532             // region of interest + topological cuts + PID
 533             Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 534
 535             //apply topological cuts
 536             if(!isTkSelected) continue;
 537             if(!isSelected) continue;
 538             if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
 539
 540             // get D candidate variables
 541             ptDStar = dstarD0pi->Pt();
 542             phiDStar = dstarD0pi->Phi();
 543             etaDStar = dstarD0pi->Eta();
 544             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
 545             if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;
 546             if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
 547             if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
 548             if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
 549
 550
 551              if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
 552
 553             phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
 554             ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
 555
 556             cout << "DStar pt = " << ptDStar << endl;
 557              cout << "pt bin = " << ptbin << endl;
 558             if(ptbin<1) continue;
 559
 560              Double_t mD0Window= fD0Window[ptbin]/3;
 561              invMassDZero = dstarD0pi->InvMassD0();
 562              deltainvMDStar = dstarD0pi->DeltaInvMass();
 563
 564
 565             // get the D meson efficiency
 566             DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
 567
 568             // check this!
 569             if(fUseDmesonEfficiencyCorrection){
 570                 if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
 571                 else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI
 572                     if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value
 573                     else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low
 574                 }
 575             }
 576             else DmesonWeight = 1.;
 577
 578
 579             // using montecarlo on reconstruction
 580             Int_t mcLabelDStar = -999;
 581             if(fmontecarlo){
 582                 // find associated MC particle for D* ->D0toKpi
 583                 mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);
 584                 if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
 585                 if(!isDStarMCtag) continue;
 586                 AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
 587                 //check if DStar from B
 588                 Int_t labelMother = MCDStar->GetMother();
 589                 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
 590                 if(!mother) continue;
 591                 Int_t motherPDG =TMath::Abs(mother->PdgCode());
 592                 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
 593             }
 594
 595
 596             // fill mass histograms
 597
 598             // cout << "Task debug check 5 " << endl;
 599             // fill D0 invariant mass
 600             if(!fmixing) {
 601                 cout << " histo name = " << Form("histDZeroMass_%d",ptbin) << endl;
 602                 ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
 603                 // cout << "Task debug check 5.1 " << endl;
 604                 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
 605             } // end if !mixing
 606
 607
 608
 609
 610
 611             // good D0 canidates
 612             if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
 613                 // fill D* invariant mass
 614                 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
 615                    // fill D* invariant mass if weighting
 616                     if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
 617                     isInPeak = kTRUE;
 618                     // fill other good candidate distributions - pt, phi eta
 619                     if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)    {
 620                         ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
 621                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 622                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 623                         nOfDStarCandidates++;
 624                         EventHasDStarCandidate=kTRUE;
 625                     } // end if in good D* mass window
 626
 627                     // count D* sideband candidates
 628                     if(fBkgMethod == kDStarSB ){
 629                         if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
 630                             ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
 631                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 632                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 633                             nOfSBCandidates++;
 634                             EventHasDStarSideBandCandidate = kTRUE;
 635                             }
 636
 637                         } // end if using DStar sidebans
 638
 639
 640             }// end good D0
 641
 642             // cout << "Task debug check 6 " << endl;
 643             // Sideband D0
 644               if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
 645                 isInDZeroSideBand = kTRUE;
 646                  if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
 647                      if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
 648
 649                      if(fBkgMethod == kDZeroSB){
 650                          if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)       {
 651
 652                              ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
 653                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 654                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 655                              nOfSBCandidates++;
 656                              EventHasDZeroSideBandCandidate = kTRUE;
 657                          }
 658                      }
 659                  }
 660             }// end SideBand D0
 661             // cout << "Task debug check 7 " << endl;
 662
 663             if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
 664             if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.151) continue; // skip all D* canidates outside the interesting mass ranges
 665            // cout << "Good D* candidate" << endl;
 666
 667             // get charge of soft pion
 668             daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
 669
 670             // get ID of daughters used for reconstruction
 671             trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
 672             trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
 673             trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
 674         }// end if reconstruction
 675
 676
 677
 678           // ============================================ using MC kinematics only
 679         Int_t DStarLabel = -1;
 680
 681         if(!fReco){ // use pure MC information
 682
 683             // get the DStar Particle
 684             DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
 685             if (!DStarMC) {
 686                 AliWarning("Careful: DStar MC Particle not found in tree, skipping");
 687                 continue;
 688             }
 689             DStarLabel = DStarMC->GetLabel();
 690             if(DStarLabel>0)isDStarMCtag = kTRUE;
 691
 692             Int_t PDG =TMath::Abs(DStarMC->PdgCode());
 693             if(PDG !=413) continue; // skip if it is not a DStar
 694
 695             // check fiducial acceptance
 696             if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
 697             ptbin=fCuts->PtBin(DStarMC->Pt());
 698
 699             //check if DStar from B
 700             Int_t labelMother = DStarMC->GetMother();
 701             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
 702             if(!mother) continue;
 703             Int_t motherPDG =TMath::Abs(mother->PdgCode());
 704             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
 705
 706
 707             Bool_t isDZero = kFALSE;
 708             Bool_t isSoftPi = kFALSE;
 709
 710             if(fUseHadronicChannelAtKineLevel){
 711                 //check decay channel on MC ************************************************
 712                 Int_t NDaugh = DStarMC->GetNDaughters();
 713                 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
 714
 715                 for(Int_t i=0; i<NDaugh;i++){ // loop on DStar daughters
 716                     Int_t daugh_label = DStarMC->GetDaughter(i);
 717                     AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
 718                     if(!mcDaughter) continue;
 719                     Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
 720                     if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
 721
 722                     if(daugh_pdg == 421) {
 723                         Int_t NDaughD0 = mcDaughter->GetNDaughters();
 724                         if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
 725                         trackiddaugh0 = mcDaughter->GetDaughter(0);
 726                         trackiddaugh1 = mcDaughter->GetDaughter(1);
 727                         Bool_t isKaon = kFALSE;
 728                         Bool_t isPion = kFALSE;
 729
 730                         for(Int_t k=0;k<NDaughD0;k++){
 731                             Int_t labelD0daugh = mcDaughter->GetDaughter(k);
 732                             AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
 733                             if(!mcGrandDaughter) continue;
 734                             Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
 735                             if(granddaugh_pdg==321) isKaon = kTRUE;
 736                             if(granddaugh_pdg==211) isPion = kTRUE;
 737                         }
 738                         if(!isKaon || !isPion) break; // skip if not correct decay channel of D0
 739                         isDZero = kTRUE;
 740                     } // end check on Dzero
 741
 742                     if(daugh_pdg == 211) {
 743                         isSoftPi = kTRUE;
 744                         daughtercharge = mcDaughter->Charge();
 745                         trackidsoftPi = daugh_label;}
 746                 } // end loop on D* daughters
 747                 if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
 748             } // end check decay channel
 749
 750             ptDStar = DStarMC->Pt();
 751             phiDStar = DStarMC->Phi();
 752             etaDStar = DStarMC->Eta();
 753
 754             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
 755
 756
 757         }// end if pure MC information
 758
 759
 760
 761
 762
 763         // getting the number of triggers in the MCtag D* case
 764         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar);
 765         if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
 766         if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
 767
 768
 769         fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
 770         fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
 771
 772
 773
 774
 775          // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
 776
 777         cout << "Correlating " << endl;
 778
 779         Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
 780
 781         // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
 782         if(fmixing && !execPool) {
 783             AliInfo("Mixed event analysis: pool is not ready");
 784             if(!isEventMixingFilledPeak && isInPeak)  {
 785                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
 786                 isEventMixingFilledPeak = kTRUE;
 787             }
 788             if (!isEventMixingFilledSB && isInDZeroSideBand)  {
 789                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
 790                 isEventMixingFilledSB=kTRUE;
 791             }
 792             continue;
 793         } // end if pool not ok
 794         // if analysis is on mixed event and pool settings are  satisfied, fill relevant histograms and continue
 795         if(fmixing&&execPool){
 796             // pool is ready - run checks on bins filling
 797             if(!isEventMixingFilledPeak && isInPeak)  {
 798                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
 799                 if(fFullmode) EventMixingChecks(aodEvent);
 800                 isEventMixingFilledPeak = kTRUE;
 801             }
 802
 803             if(!isEventMixingFilledSB && isInDZeroSideBand) {
 804                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
 805                 isEventMixingFilledSB=kTRUE;
 806             }
 807         } // end if pool ok
 808
 809
 810
 811
 812         Int_t NofEventsinPool = 1;
 813         if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
 814
 815         //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
 816
 817         for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
 818
 819             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); // process the associated tracks
 820             if(!analyzetracks) {
 821                 AliInfo("AliHFCorrelator::Cannot process the track array");
 822                 continue;
 823             }
 824
 825             Double_t arraytofill[4];
 826           // Double_t MCarraytofill[4]; // think about this
 827             Double_t weight;
 828
 829               Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
 830
 831              //************************************************** LOOP ON TRACKS *****************************************************
 832
 833              for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
 834
 835                  Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
 836                  if(!runcorrelation) continue;
 837
 838                  Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
 839                  Double_t DeltaEta = fCorrelator->GetDeltaEta();
 840
 841                  AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
 842                  if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
 843
 844                  Double_t ptHad = hadron->Pt();
 845                  Double_t phiHad = hadron->Phi();
 846                     phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
 847                  Double_t etaHad = hadron->Eta();
 848                  Int_t label = hadron->GetLabel();
 849                  Int_t trackid = hadron->GetID();
 850                  Double_t efficiency = hadron->GetWeight();
 851
 852
 853
 854
 855                  if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
 856
 857                         ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 858                         ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 859                      isTrackForPeakFilled =  kTRUE; // makes sure you do not fill twice in case of more candidates
 860                  }
 861
 862                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
 863                     ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 864                     ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 865                      isTrackForSBFilled = kTRUE;
 866                  }
 867
 868                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
 869                      ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 870                      ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 871                      isTrackForSBFilled = kTRUE;
 872                  }
 873
 874
 875                  weight = 1;
 876                  if(fUseEfficiencyCorrection && efficiency){
 877                      weight = DmesonWeight * (1./efficiency);
 878                  }
 879
 880
 881                  arraytofill[0] = DeltaPhi;
 882                  arraytofill[1] = deltainvMDStar;
 883                  arraytofill[2] = DeltaEta;
 884                  arraytofill[3] = ptHad;
 885                 // arraytofill[4] = poolbin;
 886
 887
 888                  // skip the D daughters in the correlation
 889               //   Bool_t isDdaughter = kFALSE;
 890                  if(!fmixing && fReco){ // for reconstruction
 891                      if(trackid == trackiddaugh0) continue;
 892                      if(trackid == trackiddaugh1) continue;
 893                      if(trackid == trackidsoftPi) continue;
 894                  }
 895
 896                  if(!fmixing && !fReco){ // for kinematic MC
 897                      AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
 898                      if(!part) continue;
 899                      if(IsDDaughter(DStarMC, part)) continue;
 900                      cout << "Skipping DStar  daugheter " << endl;
 901                  }
 902
 903
 904
 905                  // ============================================= FILL CORRELATION THNSparses ===============================
 906
 907                  // filling signal
 908                  if(isInPeak){
 909                      cout << "Filling signal " << endl;
 910                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 911                      //cout ("CorrelationsDStarHadron_%d",ptbin)
 912                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
 913                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
 914                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
 915                  }
 916
 917                   // filling bkg
 918                  if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
 919                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 920                       cout << "Filling bkg from D* sidebands " << endl;
 921                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
 922                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
 923                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
 924
 925                  } // end if bkg from DStar
 926
 927                  if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
 928                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 929                      cout << "Filling bkg from Dzero sidebands " << endl;
 930                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
 931                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
 932                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
 933
 934                  } // end if bkg from DZero
 935
 936
 937
 938              } // end loop on associated tracks
 939
 940         } // end loop on events in event pool
 941
 942
 943     } // end loop on D mesons
 944
 945
 946
 947
 948      Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
 949
 950
 951     if(!updated) AliInfo("Pool was not updated");
 952
 953
 954 } //end the exec
 955
 956 //________________________________________ terminate ___________________________
 957 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
 958 {
 959         // The Terminate() function is the last function to be called during
 960         // a query. It always runs on the client, it can be used to present
 961         // the results graphically or save the results to file.
 962
 963         AliAnalysisTaskSE::Terminate();
 964
 965         fOutput = dynamic_cast<TList*> (GetOutputData(1));
 966         if (!fOutput) {
 967                 printf("ERROR: fOutput not available\n");
 968                 return;
 969         }
 970
 971         return;
 972 }
 973 //_____________________________________________________
 974 Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
 975
 976     //Daughter removal in MCKine case
 977     Bool_t isDaughter = kFALSE;
 978     Int_t labelD0 = d->GetLabel();
 979
 980     Int_t mother = track->GetMother();
 981
 982     //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
 983     while (mother > 0){
 984         AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
 985         if (mcMoth){
 986             if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
 987             mother = mcMoth->GetMother(); //goes back by one
 988         } else{
 989             AliError("Failed casting the mother particle!");
 990             break;
 991         }
 992     }
 993
 994     return isDaughter;
 995 }
 996
 997 //_____________________________________________________
 998 void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
 999
1000     Double_t Pi = TMath::Pi();
1001         Int_t nbinscorr = fPhiBins;
1002         Double_t lowcorrbin = -0.5*Pi ;
1003     Double_t upcorrbin = 1.5*Pi ;
1004
1005
1006     // create ThNSparses
1007
1008     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1009
1010
1011     //sparse bins
1012
1013     //1 delta_phi
1014     //2 invariant mass D *
1015     //3 delta eta
1016     //4 track pt
1017
1018
1019     //Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
1020
1021
1022     // for reconstruction on Data and MC
1023     Int_t nbinsSparse[4]=         {nbinscorr ,     32 ,  32, 10};
1024     Double_t binLowLimitSparse[4]={lowcorrbin,0.14314 ,-1.6,  0};
1025     Double_t binUpLimitSparse[4]= {upcorrbin ,0.14794 , 1.6,  5};
1026
1027     Int_t nbinsSparseDStarSB[4]=         {nbinscorr ,     40 ,  32, 10};
1028     Double_t binLowLimitSparseDStarSB[4]={lowcorrbin,0.14788 ,-1.6,  0};
1029     Double_t binUpLimitSparseDStarSB[4]= {upcorrbin ,0.1504 , 1.6,  5};
1030
1031     TString signalSparseName = "";
1032     TString bkgSparseName = "";
1033
1034     THnSparseF * CorrelationsSignal = NULL;
1035     THnSparseF * CorrelationsBkg = NULL;
1036
1037
1038     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1039
1040
1041
1042
1043     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1044
1045
1046
1047         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1048
1049
1050
1051         signalSparseName = "CorrelationsDStar";
1052         if(fselect==1) signalSparseName += "Hadron_";
1053         if(fselect==2) signalSparseName += "Kaon_";
1054         if(fselect==3) signalSparseName += "KZero_";
1055
1056         bkgSparseName = "CorrelationsBkg";
1057         if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
1058         if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
1059         if(fselect==1) bkgSparseName += "Hadron_";
1060         if(fselect==2) bkgSparseName += "Kaon_";
1061         if(fselect==3) bkgSparseName += "KZero_";
1062
1063         signalSparseName+=Form("%d",iBin);
1064         bkgSparseName+=Form("%d",iBin);
1065         cout << "ThNSparses name = " << signalSparseName << endl;
1066
1067         // define thnsparses for signal candidates
1068         CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1069           CorrelationsSignal->Sumw2();
1070         fCorrelationOutput->Add(CorrelationsSignal);
1071
1072         // define thnsparses for bkg candidates from DStar
1073         if(fBkgMethod == kDStarSB ){
1074         CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
1075         CorrelationsBkg->Sumw2();
1076         fCorrelationOutput->Add(CorrelationsBkg);
1077         }
1078
1079         // define thnsparses for bkg candidates from DZero
1080         if(fBkgMethod == kDZeroSB ){
1081             CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1082             CorrelationsBkg->Sumw2();
1083             fCorrelationOutput->Add(CorrelationsBkg);
1084         }
1085
1086     } // end loop on bins
1087
1088
1089
1090 }
1091
1092 //__________________________________________________________________________________________________
1093 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
1094
1095     Double_t Pi = TMath::Pi();
1096         Int_t nbinscorr = fPhiBins;
1097         Double_t lowcorrbin = -0.5*Pi ;
1098     Double_t upcorrbin = 1.5*Pi ;
1099
1100     // event counter
1101     TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
1102     NofEvents->GetXaxis()->SetBinLabel(1," All events");
1103         NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1104         NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1105         NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1106         NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1107         NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1108         NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1109         NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1110     NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1111     NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1112     NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1113     fOutput->Add(NofEvents);
1114
1115     //event properties
1116     TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1117     fOutput->Add(EventPropsCheck);
1118
1119     //event properties
1120     TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",1000,0,1000);
1121     fOutput->Add(SPDMultiplicty);
1122
1123
1124     // ===================================================   D* histograms
1125     TH1F * D0mass = NULL;
1126     TH1F * DStarMass = NULL;
1127     TH1F * DStarFromSBMass = NULL;
1128
1129     TH1F * D0massWeighted = NULL;
1130     TH1F * DStarMassWeighted = NULL;
1131     TH1F * DStarFromSBMassWeighted = NULL;
1132
1133     TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "";
1134
1135     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1136     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1137
1138     //GetMinPtCandidate()
1139
1140
1141     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1142
1143
1144
1145         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1146
1147
1148         std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
1149
1150         nameDZeroMass = "histDZeroMass_";
1151         nameDStarMass = "histDStarMass_";
1152         nameDStarFromSBMass = "histDStarFromSBMass_";
1153
1154         nameDZeroMass+=Form("%d",iBin);
1155         nameDStarMass+=Form("%d",iBin);
1156         nameDStarFromSBMass+=Form("%d",iBin);
1157
1158         cout << "D zero histogram: " << nameDZeroMass << endl;
1159
1160         D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1161         DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1162         DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1163
1164         if(!fmixing){
1165         fDmesonOutput->Add(D0mass);
1166         fDmesonOutput->Add(DStarMass);
1167         fDmesonOutput->Add(DStarFromSBMass);
1168         }
1169
1170         // if using D meson efficiency, define weighted histos
1171         if(fUseDmesonEfficiencyCorrection){
1172
1173             nameDZeroMass = "histDZeroMassWeight_";
1174             nameDStarMass = "histDStarMassWeight_";
1175             nameDStarFromSBMass = "histDStarFromSBMassWeight_";
1176
1177             nameDZeroMass+=Form("%d",iBin);
1178             nameDStarMass+=Form("%d",iBin);
1179             nameDStarFromSBMass+=Form("%d",iBin);
1180
1181             D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1182             DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1183             DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1184
1185             if(!fmixing){
1186             fDmesonOutput->Add(D0massWeighted);
1187             fDmesonOutput->Add(DStarMassWeighted);
1188             fDmesonOutput->Add(DStarFromSBMassWeighted);
1189             }
1190         }
1191     }// end loop on pt bins
1192
1193
1194     // pt distributions
1195     TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
1196     fDmesonOutput->Add(RecoPtDStar);
1197         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
1198     fDmesonOutput->Add(RecoPtBkg);
1199
1200     // phi distribution
1201     TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1202     TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1203
1204     // eta distribution
1205     TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1206     TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1207
1208     if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
1209     if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
1210     if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
1211     if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
1212
1213
1214     // single track related histograms
1215     // phi distribution
1216     TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1217     TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1218
1219     // eta distribution
1220     TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1221     TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1222
1223     if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
1224     if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
1225     if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
1226     if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
1227
1228
1229     // Montecarlo for D*
1230     TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1231     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
1232
1233     TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1234     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
1235
1236         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1237         if(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
1238
1239
1240     // event mixing histograms
1241     TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1242     CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1243         CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1244     CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1245         CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1246
1247     if(fmixing) fEMOutput->Add(CheckPoolReadiness);
1248
1249
1250    /* Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
1251     Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
1252     Int_t nPoolBins = NofCentBins*NofZVrtxBins;
1253
1254
1255     TH1D * PoolBinDistribution  = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1256     fEMOutput->Add(PoolBinDistribution);
1257
1258     TH2D * EventDistributionPerPoolBin  = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1259     fEMOutput->Add(EventDistributionPerPoolBin);
1260     */
1261 }
1262
1263
1264 //__________________________________________________________________________________________________
1265 void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1266
1267
1268   //Float_t* ptbins = fCuts->GetPtBinLimits();
1269   if(fD0Window) delete fD0Window;
1270   fD0Window = new Float_t[fNofPtBins];
1271
1272   AliInfo("Enlarging the D0 mass windows from cut object\n");
1273   Int_t nvars = fCuts->GetNVars();
1274
1275   if(nvars<1){
1276     AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1277     return;
1278   }
1279   Float_t** rdcutsvalmine=new Float_t*[nvars];
1280   for(Int_t iv=0;iv<nvars;iv++){
1281     rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1282   }
1283
1284
1285     for (Int_t k=0;k<nvars;k++){
1286       for (Int_t j=0;j<fNofPtBins;j++){
1287
1288         // enlarge D0 window
1289         if(k==0)        {
1290           fD0Window[j] =fCuts->GetCutValue(0,j);
1291           rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1292           cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1293         }
1294         else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1295
1296         // set same windows
1297         //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1298       }
1299     }
1300
1301     fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1302
1303     AliInfo("\n New windows set\n");
1304     fCuts->PrintAll();
1305
1306
1307     for(Int_t iv=0;iv<nvars;iv++){
1308       delete rdcutsvalmine[iv];
1309     }
1310     delete [] rdcutsvalmine;
1311
1312 }
1313
1314
1315 //____________________________  Run checks on event mixing ___________________________________________________
1316 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
1317         // check this function
1318
1319         AliCentrality *centralityObj = 0;
1320         Int_t multiplicity = -1;
1321         Double_t MultipOrCent = -1;
1322
1323         // get the pool for event mixing
1324         if(fSystem != AA){ // pp
1325                 multiplicity = AOD->GetNTracks();
1326                 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1327         }
1328         if(fSystem == AA){ // PbPb
1329
1330                 centralityObj = AOD->GetHeader()->GetCentralityP();
1331                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1332                 AliInfo(Form("Centrality is %f", MultipOrCent));
1333         }
1334
1335         AliAODVertex *vtx = AOD->GetPrimaryVertex();
1336         Double_t zvertex = vtx->GetZ(); // zvertex
1337
1338
1339
1340
1341         AliEventPool * pool = fCorrelator->GetPool();
1342
1343
1344
1345
1346         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1347         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1348
1349         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1350         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1351 }
1352
1353
1354
1355
1356