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(7)->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     Int_t ncandidates = 0;
 489
 490     // cout << "Task debug check 3 " << endl;
 491     // loop over D meson candidates
 492     for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
 493
 494          //if(ncandidates) continue;  // if there is more than one D candidate, skip it
 495
 496         // initialize variables
 497         isInPeak = kFALSE;
 498         isInDStarSideBand = kFALSE;
 499         isInDZeroSideBand = kFALSE;
 500
 501         EventHasDStarCandidate = kFALSE;
 502         EventHasDZeroSideBandCandidate = kFALSE;
 503         EventHasDStarSideBandCandidate = kFALSE;
 504
 505
 506         isDStarMCtag = kFALSE;
 507         isDfromB = kFALSE;
 508         ptDStar = -123.4;
 509         phiDStar = -999;
 510         etaDStar = -56.;
 511         invMassDZero = - 999;
 512         deltainvMDStar = -998;
 513         AliAODRecoCascadeHF* dstarD0pi;
 514         AliAODRecoDecayHF2Prong* theD0particle;
 515         AliAODMCParticle* DStarMC=0x0;
 516         Short_t daughtercharge = -2;
 517         Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
 518         Int_t trackiddaugh1 = -1;
 519         Int_t trackidsoftPi = -1;
 520         Int_t ptbin = -1;
 521
 522         // ============================================ using reconstruction on Data or MC
 523         if(fReco){
 524             // cout << "Task debug check 4 " << endl;
 525             // get the D objects
 526             dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 527             if(!dstarD0pi->GetSecondaryVtx()) continue;
 528             theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 529             if (!theD0particle) continue;
 530
 531
 532             // apply topological cuts
 533
 534             // track quality cuts
 535             Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 536             // region of interest + topological cuts + PID
 537             Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 538
 539             //apply topological cuts
 540             if(!isTkSelected) continue;
 541             if(!isSelected) continue;
 542             if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
 543
 544             // get D candidate variables
 545             ptDStar = dstarD0pi->Pt();
 546             phiDStar = dstarD0pi->Phi();
 547             etaDStar = dstarD0pi->Eta();
 548             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
 549             if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;
 550             if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
 551             if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
 552             if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
 553
 554
 555             // if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
 556
 557             phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
 558             ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
 559
 560            // cout << "DStar pt = " << ptDStar << endl;
 561            //  cout << "pt bin = " << ptbin << endl;
 562           //  if(ptbin<3) continue;
 563
 564              Double_t mD0Window= fD0Window[ptbin]/3;
 565              invMassDZero = dstarD0pi->InvMassD0();
 566              deltainvMDStar = dstarD0pi->DeltaInvMass();
 567
 568
 569             // get the D meson efficiency
 570             DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
 571
 572             // check this!
 573             if(fUseDmesonEfficiencyCorrection){
 574                 if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
 575                 else DmesonWeight = 0; // do not consider te entry if the efficiency is too low
 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                 if(isDStarMCtag){
 587                 AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
 588                 //check if DStar from B
 589                 Int_t labelMother = MCDStar->GetMother();
 590                 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
 591                 if(!mother) continue;
 592                 Int_t motherPDG =TMath::Abs(mother->PdgCode());
 593                 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
 594                 }
 595             }
 596
 597
 598             // fill mass histograms
 599            // cout << "crash here 1" << endl;
 600             // plot D0 vs DStar mass
 601             if(!fmixing){
 602                 cout << Form("histDZerovsDStarMass_%d",ptbin) <<endl;
 603                 ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMass_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar);
 604                if(fUseDmesonEfficiencyCorrection) ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMassWeight_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar,DmesonWeight);
 605             }
 606            // cout << "crash here 2" << endl;
 607
 608
 609             // fill D0 invariant mass
 610             if(!fmixing) {
 611                 ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
 612                 // cout << "Task debug check 5.1 " << endl;
 613                 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
 614             } // end if !mixing
 615
 616
 617
 618
 619
 620             // good D0 canidates
 621             if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
 622                 // fill D* invariant mass
 623                 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
 624                    // fill D* invariant mass if weighting
 625                     if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
 626                     isInPeak = kTRUE;
 627                     // fill other good candidate distributions - pt, phi eta
 628                     if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)    {
 629                         ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
 630                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 631                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 632                         nOfDStarCandidates++;
 633                         ncandidates ++;
 634                         EventHasDStarCandidate=kTRUE;
 635                     } // end if in good D* mass window
 636
 637                     // count D* sideband candidates
 638                     if(fBkgMethod == kDStarSB ){
 639                          if(deltainvMDStar>0.1473 && deltainvMDStar<0.1644){
 640                        // if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
 641                             ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
 642                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 643                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 644                             nOfSBCandidates++;
 645                             ncandidates ++;
 646                             EventHasDStarSideBandCandidate = kTRUE;
 647                             }
 648
 649                         } // end if using DStar sidebans
 650
 651
 652             }// end good D0
 653
 654             // cout << "Task debug check 6 " << endl;
 655             // Sideband D0
 656               if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
 657                 isInDZeroSideBand = kTRUE;
 658                  if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
 659                      if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
 660
 661                      if(fBkgMethod == kDZeroSB){
 662                          if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)       {
 663
 664                              ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
 665                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
 666                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
 667                              nOfSBCandidates++;
 668                              ncandidates ++;
 669                              EventHasDZeroSideBandCandidate = kTRUE;
 670                          }
 671                      }
 672                  }
 673             }// end SideBand D0
 674             // cout << "Task debug check 7 " << endl;
 675
 676             if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
 677             if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.1644) continue; // skip all D* canidates outside the interesting mass ranges
 678            // cout << "Good D* candidate" << endl;
 679
 680             // get charge of soft pion
 681             daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
 682
 683             // get ID of daughters used for reconstruction
 684             trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
 685             trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
 686             trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
 687         }// end if reconstruction
 688
 689         // cout << "crash here 3" << endl;
 690
 691           // ============================================ using MC kinematics only
 692         Int_t DStarLabel = -1;
 693
 694         if(!fReco){ // use pure MC information
 695           //  cout << "0 - Running on kine " << endl;
 696             // get the DStar Particle
 697             DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
 698             if (!DStarMC) {
 699                 AliWarning("Careful: DStar MC Particle not found in tree, skipping");
 700                 continue;
 701             }
 702           //  cout << "1 - Have D* " << endl;
 703             DStarLabel = DStarMC->GetLabel();
 704             if(DStarLabel>0)isDStarMCtag = kTRUE;
 705
 706             Int_t PDG =TMath::Abs(DStarMC->PdgCode());
 707             if(PDG !=413) continue; // skip if it is not a DStar
 708            // cout << "PDG = " << PDG << endl;
 709             // check fiducial acceptance
 710             if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
 711           //  cout << "2 - Have D* in fiducial acceptance " << endl;
 712
 713             if(DStarMC->Pt()<fCuts->GetMinPtCandidate()||DStarMC->Pt()>fCuts->GetMaxPtCandidate()) continue;
 714             ptbin=fCuts->PtBin(DStarMC->Pt());
 715            //  cout << "3 - Have D* in ptrange acceptance " << endl;
 716             //check if DStar from B
 717             Int_t labelMother = DStarMC->GetMother();
 718             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
 719             if(!mother) continue;
 720             Int_t motherPDG =TMath::Abs(mother->PdgCode());
 721             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
 722
 723
 724             Bool_t isDZero = kFALSE;
 725             Bool_t isSoftPi = kFALSE;
 726
 727             if(fUseHadronicChannelAtKineLevel){
 728                 //check decay channel on MC ************************************************
 729                 Int_t NDaugh = DStarMC->GetNDaughters();
 730                 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
 731
 732                 for(Int_t i=0; i<NDaugh;i++){ // loop on DStar daughters
 733                     Int_t daugh_label = DStarMC->GetDaughter(i);
 734                     AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
 735                     if(!mcDaughter) continue;
 736                     Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
 737                     if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
 738
 739                     if(daugh_pdg == 421) {
 740                         Int_t NDaughD0 = mcDaughter->GetNDaughters();
 741                         if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
 742                         trackiddaugh0 = mcDaughter->GetDaughter(0);
 743                         trackiddaugh1 = mcDaughter->GetDaughter(1);
 744                         Bool_t isKaon = kFALSE;
 745                         Bool_t isPion = kFALSE;
 746
 747                         for(Int_t k=0;k<NDaughD0;k++){
 748                             Int_t labelD0daugh = mcDaughter->GetDaughter(k);
 749                             AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
 750                             if(!mcGrandDaughter) continue;
 751                             Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
 752                             if(granddaugh_pdg==321) isKaon = kTRUE;
 753                             if(granddaugh_pdg==211) isPion = kTRUE;
 754                         }
 755                         if(!isKaon || !isPion) break; // skip if not correct decay channel of D0
 756                         isDZero = kTRUE;
 757                     } // end check on Dzero
 758
 759                     if(daugh_pdg == 211) {
 760                         isSoftPi = kTRUE;
 761                         daughtercharge = mcDaughter->Charge();
 762                         trackidsoftPi = daugh_label;}
 763                 } // end loop on D* daughters
 764                 if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
 765             } // end check decay channel
 766
 767             ptDStar = DStarMC->Pt();
 768             phiDStar = DStarMC->Phi();
 769             etaDStar = DStarMC->Eta();
 770
 771             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
 772            // cout << "Dstars are selected" << endl;
 773
 774         }// end if pure MC information
 775
 776
 777
 778
 779
 780         // getting the number of triggers in the MCtag D* case
 781         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar);
 782         if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
 783         if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
 784
 785
 786         fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
 787         fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
 788
 789
 790        // cout << "crash here 4" << endl;
 791
 792          // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
 793
 794       //  cout << "Correlating " << endl;
 795
 796         Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
 797
 798         // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
 799         if(fmixing && !execPool) {
 800             AliInfo("Mixed event analysis: pool is not ready");
 801             if(!isEventMixingFilledPeak && isInPeak)  {
 802                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
 803                 isEventMixingFilledPeak = kTRUE;
 804             }
 805             if (!isEventMixingFilledSB && isInDZeroSideBand)  {
 806                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
 807                 isEventMixingFilledSB=kTRUE;
 808             }
 809             continue;
 810         } // end if pool not ok
 811         // if analysis is on mixed event and pool settings are  satisfied, fill relevant histograms and continue
 812         if(fmixing&&execPool){
 813             // pool is ready - run checks on bins filling
 814             if(!isEventMixingFilledPeak && isInPeak)  {
 815                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
 816                 if(fFullmode) EventMixingChecks(aodEvent);
 817                 isEventMixingFilledPeak = kTRUE;
 818             }
 819
 820             if(!isEventMixingFilledSB && isInDZeroSideBand) {
 821                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
 822                 isEventMixingFilledSB=kTRUE;
 823             }
 824         } // end if pool ok
 825
 826
 827
 828
 829         Int_t NofEventsinPool = 1;
 830         if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
 831
 832         Bool_t *trackOrigin = NULL;
 833        // cout << "crash here 5" << endl;
 834         //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
 835
 836         for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
 837
 838             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); // process the associated tracks
 839             if(!analyzetracks) {
 840                 AliInfo("AliHFCorrelator::Cannot process the track array");
 841                 continue;
 842             }
 843
 844             Double_t arraytofill[5];
 845             Double_t MCarraytofill[6]; // think about this
 846             Double_t weight;
 847
 848               Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
 849
 850              //************************************************** LOOP ON TRACKS *****************************************************
 851            // cout << "crash here 6" << endl;
 852              for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
 853                 // cout << "crash correlation 1" << endl;
 854                  Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
 855                  if(!runcorrelation) continue;
 856
 857                  Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
 858                  Double_t DeltaEta = fCorrelator->GetDeltaEta();
 859
 860                  AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
 861                  if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
 862
 863
 864                   Int_t trackid = hadron->GetID();
 865
 866                  // check D if it is a D meson daughter
 867                  if(!fmixing && fReco){ // for reconstruction
 868                      if(trackid == trackiddaugh0) continue;
 869                      if(trackid == trackiddaugh1) continue;
 870                      if(trackid == trackidsoftPi) continue;
 871                  }
 872
 873
 874                  Int_t label = hadron->GetLabel();
 875                  if(!fmixing && !fReco){ // for kinematic MC
 876                      AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
 877                      if(!part) continue;
 878                      if(IsDDaughter(DStarMC, part)) continue;
 879                      cout << "Skipping DStar  daugheter " << endl;
 880                  }
 881                   // if it is ok, then do the rest
 882
 883                  Double_t ptHad = hadron->Pt();
 884                  Double_t phiHad = hadron->Phi();
 885                 phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
 886                  Double_t etaHad = hadron->Eta();
 887
 888
 889                  Double_t efficiency = hadron->GetWeight();
 890
 891
 892
 893
 894                 if(fmontecarlo) trackOrigin = fAssocCuts->IsMCpartFromHF(label,fmcArray);
 895                  // cout << "crash correlation 3" << endl;
 896
 897                  if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
 898
 899                         ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 900                         ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 901                      isTrackForPeakFilled =  kTRUE; // makes sure you do not fill twice in case of more candidates
 902                  }
 903
 904                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
 905                     ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 906                     ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 907                      isTrackForSBFilled = kTRUE;
 908                  }
 909
 910                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
 911                      ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
 912                      ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
 913                      isTrackForSBFilled = kTRUE;
 914                  }
 915                  // cout << "crash correlation 4" << endl;
 916
 917                  weight = 1;
 918                  if(fUseEfficiencyCorrection && efficiency){
 919                      weight = DmesonWeight * (1./efficiency);
 920                  }
 921
 922                  // cout << "crash correlation 5" << endl;
 923                  arraytofill[0] = DeltaPhi;
 924                  arraytofill[1] = deltainvMDStar;
 925                  arraytofill[2] = DeltaEta;
 926                  arraytofill[3] = ptHad;
 927                  arraytofill[4] = poolbin;
 928
 929                  if(fmontecarlo){
 930                  MCarraytofill[0] = DeltaPhi;
 931                  MCarraytofill[1] = deltainvMDStar;
 932                  MCarraytofill[2] = DeltaEta;
 933                  MCarraytofill[3] = ptHad;
 934                  MCarraytofill[4] = poolbin;
 935
 936                 if(trackOrigin[0]) MCarraytofill[5] = 1 ;
 937                 else if(trackOrigin[1]) MCarraytofill[5] = 2 ;
 938                 else MCarraytofill[5] = 0;
 939                  }
 940
 941
 942                  // ============================================= FILL CORRELATION THNSparses ===============================
 943
 944                  // filling signal
 945                //  if(isInPeak){
 946                  if(EventHasDStarCandidate){
 947                    //  cout << "Filling signal " << endl;
 948                    // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 949                      //cout ("CorrelationsDStarHadron_%d",ptbin)
 950                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
 951                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
 952                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
 953                  }
 954
 955                   // filling bkg
 956             //     if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
 957                     if(fBkgMethod == kDStarSB && EventHasDStarSideBandCandidate){ // bkg from DStar
 958                     // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 959                     //  cout << "Filling bkg from D* sidebands " << endl;
 960                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
 961                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
 962                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
 963
 964                  } // end if bkg from DStar
 965
 966                //  if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
 967                   if(fBkgMethod == kDZeroSB && EventHasDZeroSideBandCandidate){
 968                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 969                    //  cout << "Filling bkg from Dzero sidebands " << endl;
 970                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
 971                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
 972                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
 973
 974                  } // end if bkg from DZero
 975
 976
 977
 978                  if(fmontecarlo){
 979                      // add the montecarlos
 980                        if(!fReco && TMath::Abs(etaHad)>0.8) continue;
 981                      if(!isDfromB){
 982                      //cout << "Filling correlations from charm  " << endl;
 983                      //    cout << "Ik zoek op " << Form("CorrelationsMCfromCharmHadron_%d",ptbin) << endl;
 984                      //    cout << "de lijst bevat : " << endl;
 985                          fOutputMC->ls();
 986                      if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
 987                      if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
 988                      if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
 989                      }
 990                      if(isDfromB){
 991                      //cout << "Filling correlations from beauty  " << endl;
 992                      if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
 993                      if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
 994                      if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
 995                      }
 996                  }
 997
 998
 999              } // end loop on associated tracks
1000
1001         } // end loop on events in event pool
1002
1003         if(fmixing){
1004             ((TH1D*)fEMOutput->FindObject("PoolBinDistribution"))->Fill(poolbin);
1005             ((TH2D*)fEMOutput->FindObject("EventDistributionPerPoolBin"))->Fill(poolbin,NofEventsinPool);
1006         }
1007     } // end loop on D mesons
1008
1009     if(!fmixing){
1010         if(nOfDStarCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfDCandidates"))->Fill(nOfDStarCandidates);
1011         if(nOfSBCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfSBCandidates"))->Fill(nOfSBCandidates);
1012     }
1013
1014
1015      Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
1016
1017
1018
1019
1020     if(!updated) AliInfo("Pool was not updated");
1021
1022
1023 } //end the exec
1024
1025 //________________________________________ terminate ___________________________
1026 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
1027 {
1028         // The Terminate() function is the last function to be called during
1029         // a query. It always runs on the client, it can be used to present
1030         // the results graphically or save the results to file.
1031
1032         AliAnalysisTaskSE::Terminate();
1033
1034         fOutput = dynamic_cast<TList*> (GetOutputData(1));
1035         if (!fOutput) {
1036                 printf("ERROR: fOutput not available\n");
1037                 return;
1038         }
1039
1040         return;
1041 }
1042 //_____________________________________________________
1043 Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
1044
1045     //Daughter removal in MCKine case
1046     Bool_t isDaughter = kFALSE;
1047     Int_t labelD0 = d->GetLabel();
1048
1049     Int_t mother = track->GetMother();
1050
1051     //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
1052     while (mother > 0){
1053         AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
1054         if (mcMoth){
1055             if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
1056             mother = mcMoth->GetMother(); //goes back by one
1057         } else{
1058             AliError("Failed casting the mother particle!");
1059             break;
1060         }
1061     }
1062
1063     return isDaughter;
1064 }
1065
1066 //_____________________________________________________
1067 void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
1068
1069     Double_t Pi = TMath::Pi();
1070         Int_t nbinscorr = fPhiBins;
1071         Double_t lowcorrbin = -0.5*Pi ;
1072     Double_t upcorrbin = 1.5*Pi ;
1073
1074
1075     // create ThNSparses
1076
1077     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1078
1079
1080     //sparse bins
1081
1082     //1 delta_phi
1083     //2 invariant mass D *
1084     //3 delta eta
1085     //4 track pt
1086
1087
1088     Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
1089
1090
1091     // for reconstruction on Data and MC
1092     Int_t nbinsSparse[5]=         {nbinscorr ,     32 ,  32, 10,nbinsPool};
1093     Double_t binLowLimitSparse[5]={lowcorrbin,0.14314 ,-1.6,  0,-0.5};
1094     Double_t binUpLimitSparse[5]= {upcorrbin ,0.14794 , 1.6,  5,nbinsPool-0.5};
1095
1096   //  Int_t nbinsSparseDStarSB[5]=         {nbinscorr ,     40 ,  32, 10,nbinsPool};
1097   //  Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.14788 ,-1.6,  0,-0.5};
1098   //  Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1504 , 1.6,  5,nbinsPool-0.5};
1099
1100     Int_t nbinsSparseDStarSB[5]=         {nbinscorr ,     27 ,  32, 10,nbinsPool};
1101     Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.1473 ,-1.6,  0,-0.5};
1102     Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1644 , 1.6,  5,nbinsPool-0.5};
1103
1104     Int_t nbinsSparseMC[6]=         {nbinscorr ,     32 ,  32, 10,nbinsPool,3};
1105     Double_t binLowLimitSparseMC[6]={lowcorrbin,0.14314 ,-1.6,  0,-0.5,-0.5};
1106     Double_t binUpLimitSparseMC[6]= {upcorrbin ,0.14794 , 1.6,  5,nbinsPool-0.5,2.5};
1107
1108     TString signalSparseName = "";
1109     TString bkgSparseName = "";
1110     TString MCSparseNameCharm = "";
1111     TString MCSparseNameBeauty = "";
1112
1113     THnSparseF * CorrelationsSignal = NULL;
1114     THnSparseF * CorrelationsBkg = NULL;
1115
1116     THnSparseF * MCCorrelationsCharm = NULL;
1117     THnSparseF * MCCorrelationsBeauty = NULL;
1118
1119
1120     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1121
1122
1123
1124
1125     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1126
1127
1128
1129         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1130
1131
1132
1133         signalSparseName = "CorrelationsDStar";
1134         if(fselect==1) signalSparseName += "Hadron_";
1135         if(fselect==2) signalSparseName += "Kaon_";
1136         if(fselect==3) signalSparseName += "KZero_";
1137
1138         bkgSparseName = "CorrelationsBkg";
1139         if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
1140         if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
1141         if(fselect==1) bkgSparseName += "Hadron_";
1142         if(fselect==2) bkgSparseName += "Kaon_";
1143         if(fselect==3) bkgSparseName += "KZero_";
1144
1145         MCSparseNameCharm = "CorrelationsMCfromCharm";
1146         if(fselect==1) MCSparseNameCharm += "Hadron_";
1147         if(fselect==2) MCSparseNameCharm += "Kaon_";
1148         if(fselect==3) MCSparseNameCharm += "KZero_";
1149
1150         MCSparseNameBeauty = "CorrelationsMCfromBeauty";
1151         if(fselect==1) MCSparseNameBeauty += "Hadron_";
1152         if(fselect==2) MCSparseNameBeauty += "Kaon_";
1153         if(fselect==3) MCSparseNameBeauty += "KZero_";
1154
1155         signalSparseName+=Form("%d",iBin);
1156         bkgSparseName+=Form("%d",iBin);
1157         MCSparseNameCharm+=Form("%d",iBin);
1158         MCSparseNameBeauty+=Form("%d",iBin);
1159         cout << "ThNSparses name = " << signalSparseName << endl;
1160
1161         // define thnsparses for signal candidates
1162         CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1163         CorrelationsSignal->Sumw2();
1164         fCorrelationOutput->Add(CorrelationsSignal);
1165
1166         // define thnsparses for bkg candidates from DStar
1167         if(fBkgMethod == kDStarSB ){
1168             CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",5,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
1169             CorrelationsBkg->Sumw2();
1170             fCorrelationOutput->Add(CorrelationsBkg);
1171         }
1172
1173         // define thnsparses for bkg candidates from DZero
1174         if(fBkgMethod == kDZeroSB ){
1175             CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1176             CorrelationsBkg->Sumw2();
1177             fCorrelationOutput->Add(CorrelationsBkg);
1178         }
1179
1180         if(fmontecarlo){
1181         MCCorrelationsCharm = new THnSparseF(MCSparseNameCharm.Data(),"Correlations for DStar from charm; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
1182         MCCorrelationsCharm->Sumw2();
1183         fOutputMC->Add(MCCorrelationsCharm);
1184
1185         MCCorrelationsBeauty = new THnSparseF(MCSparseNameBeauty.Data(),"Correlations for DStar from beauty; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
1186             MCCorrelationsBeauty->Sumw2();
1187             fOutputMC->Add(MCCorrelationsBeauty);
1188         }
1189
1190     } // end loop on bins
1191
1192
1193
1194 }
1195
1196 //__________________________________________________________________________________________________
1197 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
1198
1199     Double_t Pi = TMath::Pi();
1200         Int_t nbinscorr = fPhiBins;
1201         Double_t lowcorrbin = -0.5*Pi ;
1202     Double_t upcorrbin = 1.5*Pi ;
1203
1204     // event counter
1205     TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
1206     NofEvents->GetXaxis()->SetBinLabel(1," All events");
1207         NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1208         NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1209         NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1210         NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1211         NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1212         NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1213         NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1214     NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1215     NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1216     NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1217     fOutput->Add(NofEvents);
1218
1219     //event properties
1220     TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1221     fOutput->Add(EventPropsCheck);
1222
1223     //event properties
1224     TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",1000,0,1000);
1225     fOutput->Add(SPDMultiplicty);
1226
1227
1228     // ===================================================   D* histograms
1229     TH1F * D0mass = NULL;
1230     TH1F * DStarMass = NULL;
1231     TH1F * DStarFromSBMass = NULL;
1232     TH2D * DZerovsDStarMass = NULL;
1233
1234     TH1F * D0massWeighted = NULL;
1235     TH1F * DStarMassWeighted = NULL;
1236     TH1F * DStarFromSBMassWeighted = NULL;
1237     TH2D * DZerovsDStarMassWeighted = NULL;
1238
1239
1240     TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "", nameDZerovsDStarMass = "";
1241
1242     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1243     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1244
1245     //GetMinPtCandidate()
1246
1247
1248     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1249
1250
1251
1252         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1253
1254
1255         std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
1256
1257         nameDZeroMass = "histDZeroMass_";
1258         nameDStarMass = "histDStarMass_";
1259         nameDStarFromSBMass = "histDStarFromSBMass_";
1260         nameDZerovsDStarMass = "histDZerovsDStarMass_";
1261
1262         nameDZeroMass+=Form("%d",iBin);
1263         nameDStarMass+=Form("%d",iBin);
1264         nameDStarFromSBMass+=Form("%d",iBin);
1265         nameDZerovsDStarMass+=Form("%d",iBin);
1266         cout << "D vs D histogram: " << nameDZerovsDStarMass << endl;
1267
1268         D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1269         DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1270         DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1271         DZerovsDStarMass = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
1272
1273         if(!fmixing){
1274             fDmesonOutput->Add(D0mass);
1275             fDmesonOutput->Add(DStarMass);
1276             fDmesonOutput->Add(DStarFromSBMass);
1277             fDmesonOutput->Add(DZerovsDStarMass);
1278         }
1279
1280         // if using D meson efficiency, define weighted histos
1281         if(fUseDmesonEfficiencyCorrection){
1282
1283             nameDZeroMass = "histDZeroMassWeight_";
1284             nameDStarMass = "histDStarMassWeight_";
1285             nameDStarFromSBMass = "histDStarFromSBMassWeight_";
1286             nameDZerovsDStarMass = "histDZerovsDStarMassWeight_";
1287
1288             nameDZeroMass+=Form("%d",iBin);
1289             nameDStarMass+=Form("%d",iBin);
1290             nameDStarFromSBMass+=Form("%d",iBin);
1291             nameDZerovsDStarMass+=Form("%d",iBin);
1292
1293             D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1294             DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1295             DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1296             DZerovsDStarMassWeighted = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
1297
1298             if(!fmixing){
1299                 fDmesonOutput->Add(D0massWeighted);
1300                 fDmesonOutput->Add(DStarMassWeighted);
1301                 fDmesonOutput->Add(DStarFromSBMassWeighted);
1302                 fDmesonOutput->Add(DZerovsDStarMassWeighted);
1303             }
1304         }
1305     }// end loop on pt bins
1306
1307
1308     // pt distributions
1309     TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
1310     fDmesonOutput->Add(RecoPtDStar);
1311         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
1312     fDmesonOutput->Add(RecoPtBkg);
1313
1314     TH1D *NumberOfDCandidates = new TH1D("NumberOfDCandidates","Number of D candidates",10,-0.5,9.5);
1315     TH1D *NumberOfSBCandidates = new TH1D("NumberOfSBCandidates","Number of SB candidates",10,-0.5,9.5);
1316     if(!fmixing) fDmesonOutput->Add(NumberOfDCandidates);
1317     if(!fmixing) fDmesonOutput->Add(NumberOfSBCandidates);
1318
1319     // phi distribution
1320     TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1321     TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1322
1323     // eta distribution
1324     TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1325     TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1326
1327     if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
1328     if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
1329     if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
1330     if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
1331
1332
1333     // single track related histograms
1334     // phi distribution
1335     TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1336     TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1337
1338     // eta distribution
1339     TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1340     TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1341
1342     TH1D * TracksPerDcandidate = new TH1D("TracksPerDcandidate","Distribution of number of tracks per D meson candidate; N tracks; counts",200,-0.5,199.5);
1343     TH1D * TracksPerSBcandidate = new TH1D("TracksPerSBcandidate","Distribution of number of tracks per sideband candidate; N tracks; counts",200,-0.5,199.5);
1344
1345     if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
1346     if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
1347     if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
1348     if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
1349     if(!fmixing) fTracksOutput->Add(TracksPerDcandidate);
1350     if(!fmixing) fTracksOutput->Add(TracksPerSBcandidate);
1351
1352
1353     // Montecarlo for D*
1354     TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1355     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
1356
1357     TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1358     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
1359
1360         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1361         if(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
1362
1363
1364
1365     // event mixing histograms
1366     TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1367     CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1368         CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1369     CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1370         CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1371
1372     if(fmixing) fEMOutput->Add(CheckPoolReadiness);
1373
1374
1375      Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
1376      Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
1377      Int_t nPoolBins = NofCentBins*NofZVrtxBins;
1378
1379      Int_t maxevents = fAssocCuts->GetMaxNEventsInPool();
1380
1381
1382      TH1D * PoolBinDistribution  = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5);
1383    if(fmixing)  fEMOutput->Add(PoolBinDistribution);
1384
1385      TH2D * EventDistributionPerPoolBin  = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5,maxevents+2,0,maxevents+2);
1386    if(fmixing)  fEMOutput->Add(EventDistributionPerPoolBin);
1387
1388 }
1389
1390
1391 //__________________________________________________________________________________________________
1392 void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1393
1394
1395   //Float_t* ptbins = fCuts->GetPtBinLimits();
1396   if(fD0Window) delete fD0Window;
1397   fD0Window = new Float_t[fNofPtBins];
1398
1399   AliInfo("Enlarging the D0 mass windows from cut object\n");
1400   Int_t nvars = fCuts->GetNVars();
1401
1402   if(nvars<1){
1403     AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1404     return;
1405   }
1406   Float_t** rdcutsvalmine=new Float_t*[nvars];
1407   for(Int_t iv=0;iv<nvars;iv++){
1408     rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1409   }
1410
1411
1412     for (Int_t k=0;k<nvars;k++){
1413       for (Int_t j=0;j<fNofPtBins;j++){
1414
1415         // enlarge D0 window
1416         if(k==0)        {
1417           fD0Window[j] =fCuts->GetCutValue(0,j);
1418           rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1419           cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1420         }
1421         else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1422
1423         // set same windows
1424         //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1425       }
1426     }
1427
1428     fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1429
1430     AliInfo("\n New windows set\n");
1431     fCuts->PrintAll();
1432
1433
1434     for(Int_t iv=0;iv<nvars;iv++){
1435       delete rdcutsvalmine[iv];
1436     }
1437     delete [] rdcutsvalmine;
1438
1439 }
1440
1441
1442 //____________________________  Run checks on event mixing ___________________________________________________
1443 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
1444         // check this function
1445
1446         AliCentrality *centralityObj = 0;
1447         Int_t multiplicity = -1;
1448         Double_t MultipOrCent = -1;
1449
1450         // get the pool for event mixing
1451         if(fSystem != AA){ // pp
1452                 multiplicity = AOD->GetNumberOfTracks();
1453                 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1454         }
1455         if(fSystem == AA){ // PbPb
1456
1457                centralityObj = ((AliVAODHeader*)AOD->GetHeader())->GetCentralityP();
1458                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1459                 AliInfo(Form("Centrality is %f", MultipOrCent));
1460         }
1461
1462         AliAODVertex *vtx = AOD->GetPrimaryVertex();
1463         Double_t zvertex = vtx->GetZ(); // zvertex
1464
1465
1466
1467
1468         AliEventPool * pool = fCorrelator->GetPool();
1469
1470
1471
1472
1473         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1474         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1475
1476         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1477         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1478 }
1479
1480
1481
1482
1483