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