PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 //\r
  16 //\r
  17 //             Base class for DStar - Hadron Correlations Analysis\r
  18 //\r
  19 //-----------------------------------------------------------------------\r
  20 //          \r
  21 //\r
  22 //                                                 Author S.Bjelogrlic\r
  23 //                         Utrecht University \r
  24 //                      sandro.bjelogrlic@cern.ch\r
  25 //\r
  26 //-----------------------------------------------------------------------\r
  27 \r
  28 /* $Id$ */\r
  29 \r
  30 //#include <TDatabasePDG.h>\r
  31 #include <TParticle.h>\r
  32 #include <TVector3.h>\r
  33 #include <TChain.h>\r
  34 #include "TROOT.h"\r
  35 \r
  36 #include "AliAnalysisTaskDStarCorrelations.h"\r
  37 #include "AliRDHFCutsDStartoKpipi.h"\r
  38 #include "AliHFAssociatedTrackCuts.h"\r
  39 #include "AliAODRecoDecay.h"\r
  40 #include "AliAODRecoCascadeHF.h"\r
  41 #include "AliAODRecoDecayHF2Prong.h"\r
  42 #include "AliAODPidHF.h"\r
  43 #include "AliVParticle.h"\r
  44 #include "AliAnalysisManager.h"\r
  45 #include "AliAODInputHandler.h"\r
  46 #include "AliAODHandler.h"\r
  47 #include "AliESDtrack.h"\r
  48 #include "AliAODMCParticle.h"\r
  49 #include "AliNormalizationCounter.h"\r
  50 #include "AliReducedParticle.h"\r
  51 #include "AliHFCorrelator.h"\r
  52 #include "AliAODMCHeader.h"\r
  53 #include "AliEventPoolManager.h"\r
  54 #include "AliVertexingHFUtils.h"\r
  55 \r
  56 using std::cout;\r
  57 using std::endl;\r
  58 \r
  59 \r
  60 ClassImp(AliAnalysisTaskDStarCorrelations)\r
  61 \r
  62 \r
  63 //__________________________________________________________________________\r
  64 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :\r
  65 AliAnalysisTaskSE(),\r
  66 fhandler(0x0),\r
  67 fmcArray(0x0),\r
  68 fCounter(0x0),\r
  69 fCorrelator(0x0),\r
  70 fselect(0),\r
  71 fmontecarlo(kFALSE),\r
  72 fmixing(kFALSE),\r
  73 fFullmode(kFALSE),\r
  74 fSystem(pp),\r
  75 fEfficiencyVariable(kNone),\r
  76 fReco(kTRUE),\r
  77 fUseEfficiencyCorrection(kFALSE),\r
  78 fUseDmesonEfficiencyCorrection(kFALSE),\r
  79 fUseCentrality(kFALSE),\r
  80 fUseHadronicChannelAtKineLevel(kFALSE),\r
  81 fPhiBins(32),\r
  82 fEvents(0),\r
  83 fDebugLevel(0),\r
  84 fDisplacement(0),\r
  85 fDim(0),\r
  86 fNofPtBins(0),\r
  87 fMaxEtaDStar(0.9),\r
  88 fDMesonSigmas(0),\r
  89 \r
  90 fD0Window(0),\r
  91 \r
  92 fOutput(0x0),\r
  93 fOutputMC(0x0),\r
  94 fCuts(0),\r
  95 fCuts2(0),\r
  96 fUtils(0),\r
  97 fTracklets(0),\r
  98 fDeffMapvsPt(0),\r
  99 fDeffMapvsPtvsMult(0),\r
 100 fDeffMapvsPtvsEta(0)\r
 101 {\r
 102     SetDim();\r
 103     // default constructor\r
 104 }\r
 105 \r
 106 //__________________________________________________________________________\r
 107 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :\r
 108 AliAnalysisTaskSE(name),\r
 109 \r
 110 fhandler(0x0),\r
 111 fmcArray(0x0),\r
 112 fCounter(0x0),\r
 113 fCorrelator(0x0),\r
 114 fselect(0),\r
 115 fmontecarlo(kFALSE),\r
 116 fmixing(kFALSE),\r
 117 fFullmode(mode),\r
 118 fSystem(syst),\r
 119 fEfficiencyVariable(kNone),\r
 120 fReco(kTRUE),\r
 121 fUseEfficiencyCorrection(kFALSE),\r
 122 fUseDmesonEfficiencyCorrection(kFALSE),\r
 123 fUseCentrality(kFALSE),\r
 124 fUseHadronicChannelAtKineLevel(kFALSE),\r
 125 fPhiBins(32),\r
 126 fEvents(0),\r
 127 fDebugLevel(0),\r
 128 fDisplacement(0),\r
 129 fDim(0),\r
 130 fNofPtBins(0),\r
 131 fMaxEtaDStar(0.9),\r
 132 fDMesonSigmas(0),\r
 133 fD0Window(0),\r
 134 \r
 135 fOutput(0x0),\r
 136 fOutputMC(0x0),\r
 137 fCuts(0),\r
 138 fCuts2(AsscCuts),\r
 139 fUtils(0),\r
 140 fTracklets(0),\r
 141 fDeffMapvsPt(0),\r
 142 fDeffMapvsPtvsMult(0),\r
 143 fDeffMapvsPtvsEta(0)\r
 144 {\r
 145      Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");\r
 146   \r
 147     SetDim();\r
 148     if(fSystem == AA)  fUseCentrality = kTRUE; else fUseCentrality = kFALSE;\r
 149   \r
 150     fCuts=cuts;\r
 151     fNofPtBins= fCuts->GetNPtBins();\r
 152     //cout << "Enlarging the DZero window " << endl;\r
 153     EnlargeDZeroMassWindow();\r
 154    // cout << "Done" << endl;\r
 155     \r
 156    \r
 157   DefineInput(0, TChain::Class());\r
 158   \r
 159   DefineOutput(1,TList::Class()); // histos from data and MC\r
 160   DefineOutput(2,TList::Class()); // histos from MC\r
 161   DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts\r
 162   DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts\r
 163   DefineOutput(5,AliNormalizationCounter::Class());   // normalization\r
 164 }\r
 165 \r
 166 //__________________________________________________________________________\r
 167 \r
 168 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {\r
 169   //\r
 170         // destructor\r
 171         //\r
 172         \r
 173         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  \r
 174         \r
 175         if(fhandler) {delete fhandler; fhandler = 0;}\r
 176         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    \r
 177         if(fmcArray) {delete fmcArray; fmcArray = 0;}\r
 178         if(fCounter) {delete fCounter; fCounter = 0;}\r
 179         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}\r
 180         if(fOutput) {delete fOutput; fOutput = 0;}\r
 181         if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}\r
 182         if(fCuts) {delete fCuts; fCuts = 0;}\r
 183         if(fCuts2) {delete fCuts2; fCuts2=0;}\r
 184     if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}\r
 185     if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}\r
 186     if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}\r
 187 \r
 188 }\r
 189 \r
 190 //___________________________________________________________\r
 191 void AliAnalysisTaskDStarCorrelations::Init(){\r
 192         //\r
 193         // Initialization\r
 194         //\r
 195         if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");\r
 196         \r
 197         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);\r
 198         \r
 199         \r
 200     \r
 201         // Post the D* cuts\r
 202         PostData(3,copyfCuts);\r
 203         \r
 204         // Post the hadron cuts\r
 205         PostData(4,fCuts2);\r
 206     \r
 207         return;\r
 208 }\r
 209 \r
 210 \r
 211 //_________________________________________________\r
 212 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){\r
 213         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
 214         \r
 215         //slot #1\r
 216         //OpenFile(0);\r
 217         fOutput = new TList();\r
 218         fOutput->SetOwner();\r
 219         \r
 220         fOutputMC = new TList();\r
 221         fOutputMC->SetOwner();\r
 222         \r
 223         // define histograms\r
 224         DefineHistoForAnalysis();\r
 225     DefineThNSparseForAnalysis();\r
 226     \r
 227 \r
 228         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));\r
 229         fCounter->Init();\r
 230         \r
 231     Double_t Pi = TMath::Pi();\r
 232         fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb\r
 233         fCorrelator->SetDeltaPhiInterval(  -0.5*Pi, 1.5*Pi); // set correct phi interval\r
 234         //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval\r
 235         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On\r
 236         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros\r
 237         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut\r
 238         fCorrelator->SetUseMC(fmontecarlo);\r
 239         fCorrelator->SetUseReco(fReco);\r
 240     //  fCorrelator->SetKinkRemoval(kTRUE);\r
 241         Bool_t pooldef = fCorrelator->DefineEventPool();\r
 242         \r
 243         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");\r
 244     \r
 245     fUtils = new AliAnalysisUtils();\r
 246     \r
 247     \r
 248         \r
 249         PostData(1,fOutput); // set the outputs\r
 250         PostData(2,fOutputMC); // set the outputs\r
 251         PostData(5,fCounter); // set the outputs\r
 252 }\r
 253 \r
 254 //_________________________________________________\r
 255 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){\r
 256   \r
 257   \r
 258   if(fDebugLevel){\r
 259     \r
 260     if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;\r
 261     if(!fReco) std::cout << "USING MC TRUTH" << std::endl;\r
 262     std::cout << " " << std::endl;\r
 263     std::cout << "=================================================================================" << std::endl;\r
 264     if(!fmixing){\r
 265       if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;\r
 266       if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;\r
 267       if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;\r
 268     }\r
 269     if(fmixing){\r
 270       if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;\r
 271       if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;\r
 272       if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;\r
 273     }\r
 274     \r
 275   }// end if debug\r
 276   \r
 277   \r
 278   \r
 279   if (!fInputEvent) {\r
 280     Error("UserExec","NO EVENT FOUND!");\r
 281                 return;\r
 282   }\r
 283   \r
 284   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
 285   if(!aodEvent){\r
 286     AliError("AOD event not found!");\r
 287     return;\r
 288   }\r
 289   \r
 290     fTracklets = aodEvent->GetTracklets();\r
 291     \r
 292     fEvents++; // event counter\r
 293     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);\r
 294     \r
 295     fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);\r
 296     \r
 297     // load MC array\r
 298     fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
 299     if(fmontecarlo && !fmcArray){\r
 300       AliError("Array of MC particles not found");\r
 301       return;\r
 302     }\r
 303     \r
 304     \r
 305     \r
 306     \r
 307     // ********************************************** START EVENT SELECTION ****************************************************\r
 308     \r
 309     Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);\r
 310     \r
 311     if(!isEvSel) {\r
 312       \r
 313       if(fCuts->IsEventRejectedDueToPileupSPD()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);\r
 314       if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);\r
 315       if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);\r
 316       if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);\r
 317       if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);\r
 318       if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);\r
 319       if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);\r
 320       \r
 321       return;\r
 322     }\r
 323     \r
 324     // added event selection for pA\r
 325     \r
 326     if(fSystem == pA){\r
 327       \r
 328       if(fUtils->IsFirstEventInChunk(aodEvent)) {\r
 329         AliInfo("Rejecting the event - first in the chunk");\r
 330         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);\r
 331         return;\r
 332         }\r
 333       if(!fUtils->IsVertexSelected2013pA(aodEvent)) {\r
 334         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);\r
 335         AliInfo("Rejecting the event - bad vertex");\r
 336         return;\r
 337       }\r
 338     }\r
 339     // ******************************** END event selections **************************************************\r
 340     \r
 341     AliCentrality *centralityObj = 0;\r
 342     Double_t MultipOrCent = -1;\r
 343     \r
 344     if(fUseCentrality){\r
 345       /* if(fSystem == AA ){    */      centralityObj = aodEvent->GetHeader()->GetCentralityP();\r
 346       MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
 347       //AliInfo(Form("Centrality is %f", MultipOrCent));\r
 348     }\r
 349     \r
 350     if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.);\r
 351     \r
 352         \r
 353     fCorrelator->SetAODEvent(aodEvent); // set the event to be processed\r
 354     \r
 355     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);\r
 356     \r
 357     Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings\r
 358         if(!correlatorON) {\r
 359           AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");\r
 360           return;\r
 361         }\r
 362         \r
 363         if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);\r
 364                 \r
 365         // check the event type\r
 366         // load MC header\r
 367         if(fmontecarlo){\r
 368           AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));\r
 369           if (fmontecarlo && !mcHeader) {\r
 370             AliError("Could not find MC Header in AOD");\r
 371             return;\r
 372           }\r
 373           \r
 374           Bool_t isMCeventgood = kFALSE;\r
 375           \r
 376           \r
 377           Int_t eventType = mcHeader->GetEventType();\r
 378           Int_t NMCevents = fCuts2->GetNofMCEventType();\r
 379           \r
 380           for(Int_t k=0; k<NMCevents; k++){\r
 381             Int_t * MCEventType = fCuts2->GetMCEventType();\r
 382             \r
 383             if(eventType == MCEventType[k]) isMCeventgood= kTRUE;\r
 384             ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);\r
 385           }\r
 386           \r
 387           if(NMCevents && !isMCeventgood){\r
 388             if(fDebugLevel)     std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;\r
 389             return;\r
 390           }\r
 391           \r
 392         } // end if montecarlo\r
 393         \r
 394         \r
 395         // checks on vertex and multiplicity of the event\r
 396         AliAODVertex *vtx = aodEvent->GetPrimaryVertex();\r
 397         Double_t zVtxPosition = vtx->GetZ(); // zvertex\r
 398         \r
 399         \r
 400         if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);\r
 401         \r
 402         \r
 403         \r
 404         // D* reconstruction\r
 405         TClonesArray *arrayDStartoD0pi=0;\r
 406         if(!aodEvent && AODEvent() && IsStandardAOD()) {\r
 407           // In case there is an AOD handler writing a standard AOD, use the AOD\r
 408           // event in memory rather than the input (ESD) event.\r
 409           aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r
 410           // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
 411           // have to taken from the AOD event hold by the AliAODExtension\r
 412           AliAODHandler* aodHandler = (AliAODHandler*)\r
 413             ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
 414           if(aodHandler->GetExtensions()) {\r
 415             AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
 416             AliAODEvent *aodFromExt = ext->GetAOD();\r
 417             arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r
 418           }\r
 419         } else {\r
 420           arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");\r
 421         }\r
 422         \r
 423         if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
 424         \r
 425     \r
 426     // get the poolbin\r
 427     \r
 428     Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition);\r
 429   \r
 430     \r
 431         // initialize variables you will need for the D*\r
 432         \r
 433         Double_t ptDStar;//\r
 434         Double_t phiDStar;//\r
 435         Double_t etaDStar;//\r
 436         Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag;\r
 437         Double_t invMassDZero;\r
 438         Double_t deltainvMDStar;\r
 439     \r
 440         \r
 441         Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 442         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();\r
 443         \r
 444         \r
 445         //MC tagging for DStar\r
 446         //D* and D0 prongs needed to MatchToMC method\r
 447         Int_t pdgDgDStartoD0pi[2]={421,211};\r
 448         Int_t pdgDgD0toKpi[2]={321,211};\r
 449         \r
 450         Bool_t isDStarCand = kFALSE;\r
 451     Bool_t isDfromB = kFALSE;\r
 452         Bool_t isEventMixingFilledPeak = kFALSE;\r
 453         Bool_t isEventMixingFilledSB = kFALSE;\r
 454     Bool_t EventHasDStarCandidate = kFALSE;\r
 455     Bool_t EventHasDZeroSideBandCandidate = kFALSE;\r
 456     Bool_t EventHasDStarSideBandCandidate = kFALSE;\r
 457         //loop on D* candidates\r
 458         \r
 459         Int_t looponDCands = 0;\r
 460         if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast();\r
 461         if(!fReco) looponDCands = fmcArray->GetEntriesFast();\r
 462         \r
 463         Int_t nOfDStarCandidates = 0;\r
 464         Int_t nOfSBCandidates = 0;\r
 465         \r
 466         Double_t DmesonEfficiency = 1.;\r
 467         Double_t DmesonWeight = 1.;\r
 468     Double_t efficiencyvariable = -999;\r
 469     \r
 470         \r
 471         \r
 472         for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {\r
 473           isInPeak = kFALSE;\r
 474           isInDStarSideBand = kFALSE;\r
 475           isInDZeroSideBand = kFALSE;\r
 476           isDStarMCtag = kFALSE;\r
 477           isDfromB = kFALSE;\r
 478           ptDStar = -123.4;\r
 479           phiDStar = -999;\r
 480           etaDStar = -56.;\r
 481           invMassDZero = - 999;\r
 482           deltainvMDStar = -998;\r
 483           AliAODRecoCascadeHF* dstarD0pi;\r
 484           AliAODRecoDecayHF2Prong* theD0particle;\r
 485           AliAODMCParticle* DStarMC=0x0;\r
 486       Short_t daughtercharge = -2;\r
 487           Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info\r
 488           Int_t trackiddaugh1 = -1;\r
 489           Int_t trackidsoftPi = -1;\r
 490           \r
 491           // start the if reconstructed candidates from here ************************************************\r
 492           \r
 493           if(fReco){//// if reconstruction is applied\r
 494             dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);\r
 495             if(!dstarD0pi->GetSecondaryVtx()) continue;\r
 496             theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();\r
 497             if (!theD0particle) continue;\r
 498             \r
 499             \r
 500             // track quality cuts\r
 501             Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks\r
 502             // region of interest + topological cuts + PID\r
 503             Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected\r
 504             \r
 505             //apply track selections\r
 506             if(!isTkSelected) continue;\r
 507             if(!isSelected) continue;\r
 508         if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;\r
 509           \r
 510           \r
 511           ptDStar = dstarD0pi->Pt();\r
 512           phiDStar = dstarD0pi->Phi();\r
 513           etaDStar = dstarD0pi->Eta();\r
 514           if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
 515           if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;\r
 516           if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;\r
 517           if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();\r
 518           if(fEfficiencyVariable == kNone) efficiencyvariable = 0;\r
 519           \r
 520         // get the D meson efficiency\r
 521         DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);\r
 522             \r
 523          \r
 524 \r
 525             if(fUseDmesonEfficiencyCorrection){\r
 526               if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;\r
 527               else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI\r
 528                 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\r
 529                 else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low\r
 530               }\r
 531             }\r
 532             else DmesonWeight = 1.; \r
 533             \r
 534             // continue;\r
 535           \r
 536             Int_t mcLabelDStar = -999;\r
 537           if(fmontecarlo){\r
 538               // find associated MC particle for D* ->D0toKpi\r
 539               mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);\r
 540               if(mcLabelDStar>=0) isDStarMCtag = kTRUE;\r
 541               if(!isDStarMCtag) continue;\r
 542             AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);\r
 543             //check if DStar from B\r
 544             Int_t labelMother = MCDStar->GetMother();\r
 545             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
 546             if(!mother) continue;\r
 547             Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
 548             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
 549             }\r
 550                         \r
 551             \r
 552             phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar);\r
 553             \r
 554             // set the phi of the D meson in the correct range\r
 555             \r
 556             Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());\r
 557             \r
 558             \r
 559             Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma\r
 560             //  Double_t mD0Window=0.074/3;\r
 561             \r
 562             Double_t mD0Window= fD0Window[ptbin]/3;\r
 563             //cout << "Check with new window " << fD0Window[ptbin]/3 << endl;\r
 564             \r
 565             \r
 566             \r
 567             invMassDZero = dstarD0pi->InvMassD0();\r
 568             if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);\r
 569             \r
 570             deltainvMDStar = dstarD0pi->DeltaInvMass();\r
 571             \r
 572             //good D0 candidates\r
 573             if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){\r
 574               \r
 575               if(!fmixing)      ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
 576               // good D*?\r
 577               if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow){\r
 578                 \r
 579                 if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight);\r
 580                 if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);\r
 581                 isInPeak = kTRUE;\r
 582         EventHasDStarCandidate = kTRUE;\r
 583                 nOfDStarCandidates++;\r
 584               } // end Good D*\r
 585               \r
 586                 //  D* sideband?\r
 587               if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)<fDMesonSigmas[3]*dmDStarWindow)){\r
 588                 isInDStarSideBand = kTRUE;\r
 589               EventHasDStarSideBandCandidate = kTRUE;\r
 590               } // end D* sideband\r
 591               \r
 592             }// end good D0 candidates\r
 593             \r
 594             //D0 sidebands\r
 595             if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){\r
 596               if(!fmixing)((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
 597               if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight);\r
 598               \r
 599               if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0] *dmDStarWindow){ // is in DStar peak region?\r
 600                 if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight);\r
 601                 isInDZeroSideBand = kTRUE;\r
 602         EventHasDZeroSideBandCandidate = kTRUE;\r
 603                 nOfSBCandidates++;\r
 604                 if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);\r
 605               }\r
 606               \r
 607             }//end if sidebands\r
 608             \r
 609            \r
 610             \r
 611             \r
 612             if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME\r
 613             \r
 614             \r
 615             // check properties of the events containing the D*\r
 616 \r
 617      \r
 618             \r
 619             isDStarCand = kTRUE;\r
 620             \r
 621             // charge of the daughter od the\r
 622             daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();\r
 623             \r
 624             \r
 625             trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();\r
 626             trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();\r
 627             trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();\r
 628             \r
 629             // end here the reco\r
 630             \r
 631             \r
 632           }// end of if for applied reconstruction to D*\r
 633           \r
 634           Int_t DStarLabel = -1;\r
 635           \r
 636           if(!fReco){ // use pure MC information\r
 637           \r
 638         // get the DStar Particle\r
 639             DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));\r
 640             if (!DStarMC) {\r
 641               AliWarning("Careful: DStar MC Particle not found in tree, skipping");\r
 642               continue;\r
 643             }\r
 644             DStarLabel = DStarMC->GetLabel();\r
 645         if(DStarLabel>0)isDStarMCtag = kTRUE;\r
 646             \r
 647             Int_t PDG =TMath::Abs(DStarMC->PdgCode());\r
 648             if(PDG !=413) continue; // skip if it is not a DStar\r
 649             // check fiducial acceptance\r
 650         if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;\r
 651             \r
 652             //check if DStar from B\r
 653             Int_t labelMother = DStarMC->GetMother();\r
 654             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
 655                  if(!mother) continue;\r
 656             Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
 657             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
 658           \r
 659           Bool_t isDZero = kFALSE;\r
 660           Bool_t isSoftPi = kFALSE;\r
 661           \r
 662           if(fUseHadronicChannelAtKineLevel){\r
 663           //check decay channel on MC ************************************************\r
 664                 Int_t NDaugh = DStarMC->GetNDaughters();\r
 665                 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs\r
 666                 \r
 667                 for(Int_t i=0; i<NDaugh;i++){ // loop on daughters\r
 668                         Int_t daugh_label = DStarMC->GetDaughter(i);\r
 669                         AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));\r
 670                         if(!mcDaughter) continue;\r
 671                         Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());\r
 672                         if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;\r
 673               \r
 674                         if(daugh_pdg == 421) {isDZero = kTRUE;\r
 675                             Int_t NDaughD0 = mcDaughter->GetNDaughters();\r
 676                             if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs\r
 677                             trackiddaugh0 = mcDaughter->GetDaughter(0);\r
 678                             trackiddaugh1 = mcDaughter->GetDaughter(1);\r
 679                             Bool_t isKaon = kFALSE;\r
 680                             Bool_t isPion = kFALSE;\r
 681                 \r
 682                             for(Int_t k=0;k<NDaughD0;k++){\r
 683                                 Int_t labelD0daugh = mcDaughter->GetDaughter(k);\r
 684                                 AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));\r
 685                                 if(!mcGrandDaughter) continue;\r
 686                                 Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());\r
 687                                 if(granddaugh_pdg==321) isKaon = kTRUE;\r
 688                                 if(granddaugh_pdg==211) isPion = kTRUE;\r
 689                             }\r
 690                             if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0\r
 691                         }\r
 692               \r
 693                         if(daugh_pdg == 211) {\r
 694                             isSoftPi = kTRUE;\r
 695                             daughtercharge = mcDaughter->Charge();\r
 696                             trackidsoftPi = daugh_label;}\r
 697                     }\r
 698               if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel\r
 699           } // end check decay channel\r
 700             \r
 701             ptDStar = DStarMC->Pt();\r
 702             phiDStar = DStarMC->Phi();\r
 703             etaDStar = DStarMC->Eta();\r
 704           \r
 705           if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
 706             \r
 707           } // end use pure MC information\r
 708           \r
 709     \r
 710         // getting the number of triggers in the MCtag D* case\r
 711         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);\r
 712           if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);\r
 713           if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);\r
 714           \r
 715           \r
 716           fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters\r
 717           fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);\r
 718           \r
 719           \r
 720           // ************************************************ CORRELATION ANALYSIS STARTS HERE\r
 721           \r
 722           \r
 723           Bool_t execPool = fCorrelator->ProcessEventPool();\r
 724           \r
 725           \r
 726           if(fmixing && !execPool) {\r
 727             AliInfo("Mixed event analysis: pool is not ready");\r
 728             if(!isEventMixingFilledPeak && isInPeak)  {\r
 729               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1);\r
 730               isEventMixingFilledPeak = kTRUE;\r
 731             }\r
 732             if (!isEventMixingFilledSB && isInDZeroSideBand)  {\r
 733               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3);\r
 734               isEventMixingFilledSB=kTRUE;\r
 735             }\r
 736             \r
 737             continue;\r
 738           }\r
 739           \r
 740           // check event topology\r
 741           if(fmixing&&execPool){\r
 742             // pool is ready - run checks on bins filling\r
 743             if(!isEventMixingFilledPeak && isInPeak)  {\r
 744               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0);\r
 745               if(fFullmode) EventMixingChecks(aodEvent);\r
 746               isEventMixingFilledPeak = kTRUE;\r
 747             }\r
 748             \r
 749             if(!isEventMixingFilledSB && isInDZeroSideBand) {\r
 750               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2);\r
 751               isEventMixingFilledSB=kTRUE;\r
 752             }\r
 753           }\r
 754           \r
 755           Int_t NofEventsinPool = 1;\r
 756           if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();\r
 757           \r
 758           \r
 759           for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one\r
 760             \r
 761             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);\r
 762             if(!analyzetracks) {\r
 763               AliInfo("AliHFCorrelator::Cannot process the track array");\r
 764               continue;\r
 765             }\r
 766             \r
 767             //initialization of variables for correlations with leading particles\r
 768             Double_t DeltaPhiLeading = -999.;\r
 769           Double_t DeltaEtaLeading = -999.;\r
 770         \r
 771             \r
 772             Int_t NofTracks = fCorrelator->GetNofTracks();\r
 773             \r
 774             \r
 775             if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks);\r
 776             if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks);\r
 777             \r
 778             \r
 779             \r
 780             Double_t arraytofill[5];\r
 781             Double_t MCarraytofill[7];\r
 782             \r
 783             \r
 784             Double_t weight;\r
 785             \r
 786           for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates\r
 787               Bool_t runcorrelation = fCorrelator->Correlate(iTrack);\r
 788               if(!runcorrelation) continue;\r
 789               \r
 790               Double_t DeltaPhi = fCorrelator->GetDeltaPhi();\r
 791               Double_t DeltaEta = fCorrelator->GetDeltaEta();\r
 792               \r
 793               AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();\r
 794               if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}\r
 795               \r
 796               Double_t ptHad = hadron->Pt();\r
 797               Double_t phiHad = hadron->Phi();\r
 798               Double_t etaHad = hadron->Eta();\r
 799               Int_t label = hadron->GetLabel();\r
 800               Int_t trackid = hadron->GetID();\r
 801               Double_t efficiency = hadron->GetWeight();\r
 802               \r
 803               weight = 1;\r
 804               if(fUseEfficiencyCorrection && efficiency){\r
 805                   weight = DmesonWeight * (1./efficiency);\r
 806               }\r
 807         \r
 808               phiHad = fCorrelator->SetCorrectPhiRange(phiHad);\r
 809               \r
 810               \r
 811               if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency);\r
 812               \r
 813               arraytofill[0] = DeltaPhi;\r
 814               arraytofill[1] = DeltaEta;\r
 815               arraytofill[2] = ptDStar;\r
 816               arraytofill[3] = ptHad;\r
 817               arraytofill[4] = poolbin;\r
 818                 \r
 819               \r
 820               MCarraytofill[0] = DeltaPhi;\r
 821               MCarraytofill[1] = DeltaEta;\r
 822               MCarraytofill[2] = ptDStar;\r
 823               MCarraytofill[3] = ptHad;\r
 824               MCarraytofill[4] = poolbin;\r
 825               \r
 826              \r
 827               if(fmontecarlo){\r
 828                 if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks);\r
 829                 if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks);\r
 830                 if(label<0) continue; // skip track with wrong label\r
 831               }\r
 832               \r
 833                Bool_t isDdaughter = kFALSE;\r
 834             // skip the D daughters in the correlation\r
 835               if(!fmixing && fReco){\r
 836                   if(trackid == trackiddaugh0) continue;\r
 837                   if(trackid == trackiddaugh1) continue;\r
 838                   if(trackid == trackidsoftPi) continue;\r
 839               }\r
 840               \r
 841               if(!fmixing && !fReco){\r
 842                   AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);\r
 843                   if(!part) continue;\r
 844                   if(IsDDaughter(DStarMC, part)) continue;\r
 845                   cout << "Skipping DStar  daugheter " << endl;\r
 846               }\r
 847               if(!fmixing && !fReco && fmontecarlo){  // skip D* Daughetrs if it is Pure MCDStar\r
 848                     Int_t hadronlabel = label;\r
 849                     for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers\r
 850                         if(DStarLabel<0){ break;}\r
 851                         if(hadronlabel<0) { break;}\r
 852                         AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fmcArray->At(hadronlabel));\r
 853                         if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;}\r
 854                         hadronlabel = mcParticle->GetMother();\r
 855                         if(hadronlabel == DStarLabel) isDdaughter = kTRUE;\r
 856                     }\r
 857                 \r
 858                   if(isDdaughter && fDebugLevel){\r
 859                       std::cout << "It is the D* daughter with label " << label << std::endl;\r
 860                       std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl;\r
 861                       std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl;\r
 862                       std::cout << "Soft pi label = " << trackidsoftPi << std::endl;\r
 863                   }\r
 864                     \r
 865                                         if(isDdaughter) continue; // skip if track is from DStar\r
 866                                 }\r
 867                 \r
 868                                 // ================ FILL CORRELATION HISTOGRAMS ===============================\r
 869                                 \r
 870                 // monte carlo case (mc tagged D*)\r
 871               if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo\r
 872                 \r
 873                 Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)\r
 874                               \r
 875                 MCarraytofill[5] = 0;\r
 876                 if(PartSource[0]) MCarraytofill[5] = 1;\r
 877                 if(PartSource[1]) MCarraytofill[5] = 2;\r
 878                 if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3;\r
 879                     if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4;\r
 880                     if(PartSource[3]) MCarraytofill[5] = 5;\r
 881                     if(!isDfromB) MCarraytofill[6] = 0;\r
 882                     if(isDfromB) MCarraytofill[6] = 1;\r
 883                     if(!fReco && TMath::Abs(etaHad)>0.8) {\r
 884                       delete [] PartSource;\r
 885                       continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
 886                     }\r
 887                     ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill);\r
 888                     \r
 889                     delete[] PartSource;\r
 890                                 }\r
 891               \r
 892                 // Good DStar canidates\r
 893                                 if(isInPeak)  {\r
 894                     \r
 895                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
 896                     if(fselect==1)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
 897                     if(fselect==2)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
 898                     if(fselect==3)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
 899                                         \r
 900                                     ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent);\r
 901                                         if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad);\r
 902                                     \r
 903                                 }\r
 904               \r
 905                 // Sidebands from D0 candidate\r
 906                                 if(isInDZeroSideBand) {\r
 907                                         \r
 908                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
 909                     if(fselect==1)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
 910                     if(fselect==2)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
 911                     if(fselect==3)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
 912                     \r
 913                                         if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad);\r
 914                     \r
 915                 }\r
 916               \r
 917               // Sidebands from D* candidate\r
 918                 if(isInDStarSideBand) {\r
 919                                         \r
 920                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
 921                     if(fselect==1 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
 922                     if(fselect==2 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
 923                     if(fselect==3 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
 924 \r
 925                                 }\r
 926                 \r
 927                 \r
 928                         } // end loop on track candidates\r
 929             \r
 930                         \r
 931                         \r
 932             // fill the leading particle histograms\r
 933             \r
 934             if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
 935                         if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
 936             \r
 937                 } // end loop on events in the pool\r
 938         \r
 939         }// end loop on D* candidates\r
 940         \r
 941     \r
 942  \r
 943         // check events with D* or SB canidates\r
 944     if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition);\r
 945      if(fFullmode && EventHasDZeroSideBandCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition);\r
 946     \r
 947     if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition);\r
 948     \r
 949     \r
 950         if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent);\r
 951     if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent);\r
 952         \r
 953     // update event pool\r
 954     Bool_t updated = fCorrelator->PoolUpdate();\r
 955         \r
 956     //  if(updated) EventMixingChecks(aodEvent);\r
 957         if(!updated) AliInfo("Pool was not updated");\r
 958         \r
 959     \r
 960 } //end the exec\r
 961 \r
 962 //________________________________________ terminate ___________________________\r
 963 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)\r
 964 {    \r
 965         // The Terminate() function is the last function to be called during\r
 966         // a query. It always runs on the client, it can be used to present\r
 967         // the results graphically or save the results to file.\r
 968         \r
 969         AliAnalysisTaskSE::Terminate();\r
 970         \r
 971         fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
 972         if (!fOutput) {     \r
 973                 printf("ERROR: fOutput not available\n");\r
 974                 return;\r
 975         }\r
 976 \r
 977         return;\r
 978 }\r
 979 //_____________________________________________________\r
 980 Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {\r
 981     \r
 982     //Daughter removal in MCKine case\r
 983     Bool_t isDaughter = kFALSE;\r
 984     Int_t labelD0 = d->GetLabel();\r
 985     \r
 986     Int_t mother = track->GetMother();\r
 987     \r
 988     //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)\r
 989     while (mother > 0){\r
 990         AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!\r
 991         if (mcMoth){\r
 992             if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;\r
 993             mother = mcMoth->GetMother(); //goes back by one\r
 994         } else{\r
 995             AliError("Failed casting the mother particle!");\r
 996             break;\r
 997         }\r
 998     }\r
 999     \r
1000     return isDaughter;\r
1001 }\r
1002 \r
1003 //_____________________________________________________\r
1004 void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){\r
1005     \r
1006     //cout << "DEFINING THNSPARSES "<< endl;\r
1007     \r
1008     Double_t Pi = TMath::Pi();\r
1009         Int_t nbinscorr = fPhiBins;\r
1010         Double_t lowcorrbin = -0.5*Pi;\r
1011         Double_t upcorrbin = 1.5*Pi;\r
1012     // define the THnSparseF\r
1013     \r
1014     //sparse bins\r
1015     \r
1016     //1 delta_phi\r
1017     //2 delta_eta\r
1018     //3 D* pt\r
1019     //4 multiplicity\r
1020     //5 track pt\r
1021     //6 zVtx position\r
1022     \r
1023     Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins());\r
1024     \r
1025     \r
1026     Int_t nbinsSparse[5]=         {nbinscorr,   32,30,250,nbinsPool};\r
1027     Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0,  0,-0.5};\r
1028     Double_t binUpLimitSparse[5]= {upcorrbin,  1.6,30,25,nbinsPool-0.5};\r
1029   \r
1030     Int_t MCnbinsSparse[7]=         {nbinscorr,   32,30,250,nbinsPool,10,2};    \r
1031     Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0,  0,-0.5,-0.5,-0.5};      //  \r
1032     Double_t MCbinUpLimitSparse[7]= {upcorrbin,  1.6,30,25,nbinsPool-0.5,9.5,1.5};\r
1033     \r
1034     TString sparsename = "CorrelationsDStar";\r
1035     if(fselect==1) sparsename += "Hadron";\r
1036         if(fselect==2) sparsename += "Kaon";\r
1037         if(fselect==3) sparsename += "KZero";\r
1038     \r
1039     TString D0Bkgsparsename = "DZeroBkg";\r
1040     D0Bkgsparsename += sparsename;\r
1041     \r
1042     TString DStarBkgsparsename = "DStarBkg";\r
1043     DStarBkgsparsename += sparsename;\r
1044     \r
1045     TString MCSparseName = "MCDStar";\r
1046     MCSparseName += sparsename;\r
1047     // signal correlations\r
1048     THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
1049     \r
1050     // bkg correlations from D0 sidebands\r
1051     THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
1052     \r
1053     // bkg correlations from D* sidebands\r
1054     THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
1055     \r
1056     \r
1057     THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse);\r
1058     \r
1059     MCCorrelations->GetAxis(5)->SetBinLabel(1," All ");\r
1060         MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour");\r
1061         MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D");\r
1062         MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D");\r
1063         MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B");\r
1064         MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour");\r
1065         MCCorrelations->GetAxis(5)->SetBinLabel(7," from c");\r
1066         MCCorrelations->GetAxis(5)->SetBinLabel(8," from b");\r
1067     \r
1068     MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c");\r
1069     MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b");\r
1070     \r
1071     Correlations->Sumw2();\r
1072     DZeroBkgCorrelations->Sumw2();\r
1073     DStarBkgCorrelations->Sumw2();\r
1074     \r
1075     fOutput->Add(Correlations);\r
1076     fOutput->Add(DZeroBkgCorrelations);\r
1077     if(fFullmode) fOutput->Add(DStarBkgCorrelations);\r
1078     if(fmontecarlo) fOutputMC->Add(MCCorrelations);\r
1079     \r
1080     \r
1081 }\r
1082 //__________________________________________________________________________________________________\r
1083 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){\r
1084         \r
1085         Double_t Pi = TMath::Pi();\r
1086         Int_t nbinscorr = fPhiBins;\r
1087         Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center\r
1088         Double_t upcorrbin = 1.5*Pi ;\r
1089         \r
1090         // ========================= histograms for both Data and MonteCarlo\r
1091         \r
1092         \r
1093         TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);\r
1094     NofEvents->GetXaxis()->SetBinLabel(1," All events");\r
1095         NofEvents->GetXaxis()->SetBinLabel(2," Selected events");\r
1096         NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");\r
1097         NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");\r
1098         NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");\r
1099         NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");\r
1100         NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");\r
1101         NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");\r
1102     NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");\r
1103     NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");\r
1104     NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");\r
1105     fOutput->Add(NofEvents);\r
1106         \r
1107         \r
1108         \r
1109         \r
1110         TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);\r
1111         if(!fmixing && fFullmode) fOutput->Add(D0InvMass);\r
1112         \r
1113         TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);\r
1114         if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB);\r
1115         \r
1116         //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);\r
1117         //if(!fmixing) fOutput->Add(DeltaInvMass);\r
1118     TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution; D* p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);\r
1119         if(!fmixing) fOutput->Add(DeltaInvMass);\r
1120         \r
1121         TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution; SB p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);\r
1122         if(!fmixing) fOutput->Add(bkgDeltaInvMass);\r
1123     \r
1124     DeltaInvMass->Sumw2();\r
1125     bkgDeltaInvMass->Sumw2();\r
1126         \r
1127         TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);\r
1128         if(!fmixing) fOutput->Add(RecoPtDStar);\r
1129         \r
1130         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);\r
1131         if(!fmixing) fOutput->Add(RecoPtBkg);\r
1132     \r
1133     TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);\r
1134     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);\r
1135     \r
1136     TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);\r
1137     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);\r
1138         \r
1139         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);\r
1140         if(!fmixing) fOutput->Add(MCtagPtDStar);\r
1141         \r
1142         TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);\r
1143     if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra);\r
1144         \r
1145         TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);\r
1146         if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF);\r
1147         \r
1148         TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5);\r
1149         if(fFullmode) fOutput->Add(NofTracksInPeak);\r
1150         \r
1151         TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5);\r
1152         if(fFullmode) fOutput->Add(NofTracksInSB);\r
1153         \r
1154         TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
1155         if(fFullmode)fOutput->Add(TracksInPeakSpectra);\r
1156         \r
1157         TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
1158         if(fFullmode)fOutput->Add(TracksInSBSpectra);\r
1159         \r
1160         \r
1161         //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);\r
1162         //if(fmixing) fOutput->Add(EventMixingCheck);\r
1163     \r
1164     \r
1165     TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
1166         if(fFullmode)fOutput->Add(EventPropsCheck);\r
1167     \r
1168     TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
1169         if(fFullmode)fOutput->Add(EventPropsCheckifDStar);\r
1170     \r
1171     TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
1172         if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB);\r
1173     \r
1174     TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
1175         if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB);\r
1176     \r
1177         \r
1178     TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005);\r
1179         if(fFullmode)fOutput->Add(WeightChecks);\r
1180     \r
1181         \r
1182         \r
1183         TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
1184         fOutput->Add(PhiEtaTrigger);\r
1185         \r
1186         TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
1187         fOutput->Add(PhiEtaSideBand);\r
1188         \r
1189         TH3F * PhiEtaPart = new TH3F("PhiEtaPart","#phi distribution of the associated particle; #phi; #eta; multiplicity",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9,100,0,1000);\r
1190         fOutput->Add(PhiEtaPart);\r
1191     \r
1192     TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
1193         if(fFullmode)fOutput->Add(DStarCandidates);\r
1194     \r
1195     TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
1196         if(fFullmode)fOutput->Add(SBCandidates);\r
1197         \r
1198         \r
1199         //correlations histograms\r
1200         TString histoname1 = "DPhiDStar";\r
1201         if(fselect==1) histoname1 += "Hadron";\r
1202         if(fselect==2) histoname1 += "Kaon";\r
1203         if(fselect==3) histoname1 += "KZero";\r
1204         \r
1205         /*\r
1206      TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1207      \r
1208      TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1209      \r
1210      //side band background histograms\r
1211      TString histoname2 = "bkg";\r
1212      histoname2 += histoname1;\r
1213      TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1214      TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1215      \r
1216      \r
1217      fOutput->Add(DPhiDStar);\r
1218      \r
1219      if(fselect==3){fOutput->Add(DPhiDStarKZero1);}\r
1220      \r
1221      fOutput->Add(bkgDPhiDStar);\r
1222      \r
1223      if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}\r
1224      \r
1225      */\r
1226         // leading particle\r
1227         TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1228         TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1229         \r
1230         if(fFullmode)fOutput->Add(leadingcand);\r
1231         if(fFullmode)fOutput->Add(leadingsidebands);\r
1232         \r
1233         // ========================= histos for analysis on MC only\r
1234         \r
1235         TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);\r
1236         if(fmontecarlo) fOutputMC->Add(EventTypeMC);\r
1237         \r
1238         TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);\r
1239         MCSources->GetXaxis()->SetBinLabel(1," All ");\r
1240         MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
1241         MCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
1242         MCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
1243         MCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
1244         MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
1245         MCSources->GetXaxis()->SetBinLabel(7," from c");\r
1246         MCSources->GetXaxis()->SetBinLabel(8," from b");\r
1247         \r
1248         if(fmontecarlo) fOutputMC->Add(MCSources);\r
1249     \r
1250     // leading particle from mc source\r
1251     TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);\r
1252         LeadingMCSources->GetXaxis()->SetBinLabel(1," All ");\r
1253         LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
1254         LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
1255         LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
1256         LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
1257         LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
1258         LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");\r
1259         LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");\r
1260         \r
1261         if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources);\r
1262         \r
1263     // all hadrons\r
1264         TString histoname3 = "MCTag";\r
1265         histoname3 += histoname1;\r
1266         TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1267         \r
1268         TString histoname44 = "CharmDOrigin";\r
1269         histoname44 += histoname1;\r
1270         histoname44 += "MC";\r
1271         \r
1272         TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1273         \r
1274         \r
1275         TString histoname54 = "BeautyDOrigin";\r
1276         histoname54 += histoname1;\r
1277         histoname54 += "MC";\r
1278         TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1279         \r
1280         TString histoname55 = "BeautyBOrigin";\r
1281         histoname55 += histoname1;\r
1282         histoname55 += "MC";\r
1283         TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1284         \r
1285         TString histoname4 = "CharmQuarkOrigin";\r
1286         histoname4 += histoname1;\r
1287         histoname4 += "MC";\r
1288         TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1289         \r
1290         TString histoname5 = "BeautyQuarkOrigin";\r
1291         histoname5 += histoname1;\r
1292         histoname5 += "MC";\r
1293         TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1294         \r
1295         TString histoname6 = "NonHFOrigin";\r
1296         histoname6 += histoname1;\r
1297         histoname6 += "MC";\r
1298         TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1299     \r
1300     if(fmontecarlo && fFullmode){\r
1301         \r
1302         fOutputMC->Add(MCTagDPhiDStar);\r
1303         fOutputMC->Add(CharmDOriginDPhiDStar);\r
1304         fOutputMC->Add(BeautyDOriginDPhiDStar);\r
1305         fOutputMC->Add(BeautyBOriginDPhiDStar);\r
1306         fOutputMC->Add(CharmQuarkOriginDPhiDStar);\r
1307         fOutputMC->Add(BeautyQuarkOriginDPhiDStar);\r
1308                 fOutputMC->Add(NonHFOriginDPhiDStar);\r
1309         \r
1310         }\r
1311     \r
1312     // ========================= histos for analysis on MC\r
1313     // all leading hadron\r
1314         TString Leadinghistoname3 = "LeadingMCTag";\r
1315         Leadinghistoname3 += histoname1;\r
1316         TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1317     \r
1318         TString Leadinghistoname44 = "LeadingCharmDOrigin";\r
1319         Leadinghistoname44 += histoname1;\r
1320         Leadinghistoname44 += "MC";\r
1321         \r
1322         TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1323         \r
1324         \r
1325         TString Leadinghistoname54 = "LeadingBeautyDOrigin";\r
1326         Leadinghistoname54 += histoname1;\r
1327         Leadinghistoname54 += "MC";\r
1328         TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1329         \r
1330         TString Leadinghistoname55 = "LeadingBeautyBOrigin";\r
1331         Leadinghistoname55 += histoname1;\r
1332         Leadinghistoname55 += "MC";\r
1333         TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1334         \r
1335         TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";\r
1336         Leadinghistoname4 += histoname1;\r
1337         Leadinghistoname4 += "MC";\r
1338         TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1339         \r
1340         TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";\r
1341         Leadinghistoname5 += histoname1;\r
1342         Leadinghistoname5 += "MC";\r
1343         TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
1344     \r
1345     \r
1346         \r
1347         \r
1348         if(fmontecarlo && fFullmode){\r
1349                 \r
1350                 fOutputMC->Add(LeadingMCTagDPhiDStar);\r
1351                 fOutputMC->Add(LeadingCharmDOriginDPhiDStar);\r
1352                 fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);\r
1353                 fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);\r
1354                 fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);\r
1355                 fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);\r
1356                 \r
1357         }\r
1358         \r
1359         TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5);\r
1360         MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");\r
1361         MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");\r
1362         MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");\r
1363         MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");\r
1364         MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");\r
1365         MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");\r
1366         if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);\r
1367         \r
1368         TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5);\r
1369         if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels);\r
1370         \r
1371         // ============================= EVENT MIXING CHECKS ======================================\r
1372         \r
1373         Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();\r
1374         Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();\r
1375         Int_t NofCentBins = fCuts2->GetNCentPoolBins();\r
1376         Double_t * CentBins = fCuts2->GetCentPoolBins();\r
1377         Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();\r
1378         Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();\r
1379     \r
1380     \r
1381     \r
1382         Int_t k =0;\r
1383         \r
1384         if(fSystem == AA) k = 100; // PbPb centrality\r
1385         if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity\r
1386         \r
1387         \r
1388         //Double_t minvalue = CentBins[0];\r
1389         //Double_t maxvalue = CentBins[NofCentBins+1];\r
1390         //Double_t Zminvalue = ZVrtxBins[0];\r
1391         //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];\r
1392         \r
1393         Double_t minvalue, maxvalue;\r
1394     Double_t Zminvalue, Zmaxvalue;\r
1395     \r
1396     Zminvalue = -15.;\r
1397     Zmaxvalue = 15;\r
1398     if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb\r
1399     if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity\r
1400     \r
1401         //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
1402     // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
1403     //  Double_t * events = Nevents;\r
1404     Double_t eventsv[] ={0,1000000};\r
1405     //Double_t * events = new Double_t[2];\r
1406    // events[0] = 0;\r
1407 //      events[1] = 1000000;\r
1408     Double_t *events = eventsv;\r
1409     Int_t Nevents = 1000000;\r
1410     //  TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events);\r
1411     \r
1412     TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]);\r
1413     \r
1414         EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
1415         EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
1416         EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");\r
1417         if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin);\r
1418         \r
1419         Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;\r
1420         //Int_t Diff = MaxNofTracks-MinNofTracks;\r
1421         \r
1422     //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};\r
1423     //  Double_t  * trackN = Ntracks;\r
1424         \r
1425         TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks);\r
1426         NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
1427         NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
1428         NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");\r
1429         \r
1430         if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin);\r
1431         \r
1432         TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);\r
1433         NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
1434         NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");\r
1435         if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls);\r
1436         \r
1437     \r
1438         \r
1439         TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);\r
1440         EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
1441         EventProps->GetYaxis()->SetTitle("Z vertex [cm]");\r
1442         if(fmixing && fFullmode) fOutput->Add(EventProps);\r
1443     \r
1444     TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);\r
1445     CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");\r
1446         CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");\r
1447     CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");\r
1448         CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");\r
1449         \r
1450     if(fmixing) fOutput->Add(CheckPoolReadiness);\r
1451     \r
1452         \r
1453 }\r
1454 \r
1455 //__________________________________________________________________________________________________\r
1456 void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){\r
1457     \r
1458 \r
1459   //Float_t* ptbins = fCuts->GetPtBinLimits();\r
1460   if(fD0Window) delete fD0Window;\r
1461   fD0Window = new Float_t[fNofPtBins];\r
1462     \r
1463   AliInfo("Enlarging the D0 mass windows from cut object\n"); \r
1464   Int_t nvars = fCuts->GetNVars();\r
1465 \r
1466   if(nvars<1){\r
1467     AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");\r
1468     return;\r
1469   }\r
1470   Float_t** rdcutsvalmine=new Float_t*[nvars];\r
1471   for(Int_t iv=0;iv<nvars;iv++){\r
1472     rdcutsvalmine[iv]=new Float_t[fNofPtBins];\r
1473   }\r
1474     \r
1475     \r
1476     for (Int_t k=0;k<nvars;k++){\r
1477       for (Int_t j=0;j<fNofPtBins;j++){\r
1478             \r
1479         // enlarge D0 window\r
1480         if(k==0)        {\r
1481           fD0Window[j] =fCuts->GetCutValue(0,j);\r
1482           rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);\r
1483           cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;\r
1484         }\r
1485         else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);\r
1486                         \r
1487         // set same windows\r
1488         //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);\r
1489       }\r
1490     }\r
1491     \r
1492     fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);\r
1493     \r
1494     AliInfo("\n New windows set\n");     \r
1495     fCuts->PrintAll();\r
1496     \r
1497     \r
1498     for(Int_t iv=0;iv<nvars;iv++){\r
1499       delete rdcutsvalmine[iv];\r
1500     }\r
1501     delete [] rdcutsvalmine;\r
1502     \r
1503 }\r
1504 \r
1505 \r
1506 //____________________________  Run checks on event mixing ___________________________________________________\r
1507 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){\r
1508         \r
1509     \r
1510         AliCentrality *centralityObj = 0;\r
1511         Int_t multiplicity = -1;\r
1512         Double_t MultipOrCent = -1;\r
1513         \r
1514         // get the pool for event mixing\r
1515         if(fSystem != AA){ // pp\r
1516                 multiplicity = AOD->GetNTracks();\r
1517                 MultipOrCent = multiplicity; // convert from Int_t to Double_t\r
1518         }\r
1519         if(fSystem == AA){ // PbPb\r
1520                 \r
1521                 centralityObj = AOD->GetHeader()->GetCentralityP();\r
1522                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
1523                 AliInfo(Form("Centrality is %f", MultipOrCent));\r
1524         }\r
1525         \r
1526         AliAODVertex *vtx = AOD->GetPrimaryVertex();\r
1527         Double_t zvertex = vtx->GetZ(); // zvertex\r
1528         \r
1529         \r
1530         \r
1531         \r
1532         AliEventPool * pool = fCorrelator->GetPool();\r
1533         \r
1534     \r
1535         \r
1536         \r
1537         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool\r
1538         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties\r
1539         \r
1540         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool\r
1541         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool\r
1542 }\r
1543         \r
1544 \r
1545 \r
1546 \r
1547 \r