PWGHF/correlationHF/AliAnalysisTaskSEDplusCorrelations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*************************************************************************
  19 // Class AliAnalysisTaskSEDplusCorrelations
  20 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
  21 //comparison of heavy-flavour decay candidates
  22 // to MC truth (kinematics stored in the AOD)
  23 // Authors: Sadhana Dash (correlation)
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <TClonesArray.h>
  27 #include <TNtuple.h>
  28 #include <TCanvas.h>
  29 #include <TList.h>
  30 #include <TString.h>
  31 #include <TDatabasePDG.h>
  32
  33 #include <AliAnalysisDataSlot.h>
  34 #include <AliAnalysisDataContainer.h>
  35 #include "AliAnalysisManager.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODPidHF.h"
  40 #include "AliAODVertex.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODRecoDecayHF3Prong.h"
  43 #include "AliAnalysisVertexingHF.h"
  44 #include "AliAnalysisTaskSE.h"
  45 #include "AliAnalysisTaskSEDplusCorrelations.h"
  46 #include "AliNormalizationCounter.h"
  47 #include "AliVParticle.h"
  48 #include "AliHFAssociatedTrackCuts.h"
  49 #include "AliReducedParticle.h"
  50 #include "AliHFCorrelator.h"
  51
  52
  53 using std::cout;
  54 using std::endl;
  55
  56 ClassImp(AliAnalysisTaskSEDplusCorrelations)
  57
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskSEDplusCorrelations::AliAnalysisTaskSEDplusCorrelations():
  61 AliAnalysisTaskSE(),
  62   fOutput(0),
  63   fCorrelator(0),
  64   fSelect(0),
  65   fDisplacement(0),
  66   fHistNEvents(0),
  67   fMCSources(0),
  68   fK0Origin(0),
  69   fKaonOrigin(0),
  70   fInvMassK0S(0),
  71   fEventMix(0),
  72   fPtVsMass(0),
  73   fPtVsMassTC(0),
  74   fYVsPt(0),
  75   fYVsPtTC(0),
  76   fYVsPtSig(0),
  77   fYVsPtSigTC(0),
  78   fUpmasslimit(1.965),
  79   fLowmasslimit(1.765),
  80   fNPtBins(0),
  81   fBinWidth(0.002),
  82   fListCuts(0),
  83   fListCutsAsso(0),
  84   fRDCutsAnalysis(0),
  85   fCuts(0),
  86   fCounter(0),
  87   fReadMC(kFALSE),
  88   fUseBit(kTRUE),
  89   fMixing(kFALSE),
  90   fSystem(kFALSE)
  91 {
  92   // Default constructor
  93
  94   for(Int_t i=0;i<3*kMaxPtBins;i++){
  95     fMassVsdPhiHistHad[i]=0;
  96     fMassVsdEtaHistHad[i]=0;
  97     fMassVsdPhiHistKaon[i]=0;
  98     fMassVsdEtaHistKaon[i]=0;
  99     fMassVsdPhiHistKshort[i]=0;
 100     fMassVsdEtaHistKshort[i]=0;
 101     fMassVsdPhiHistLeadHad[i]=0;
 102     fMassVsdEtaHistLeadHad[i]=0;
 103     fMassHistK0S[i]=0;
 104     fLeadPt[i]=0;
 105     fMassHist[i]=0;
 106     fMassHistTC[i]=0;
 107     fMassHistTCPlus[i]=0;
 108     fMassHistTCMinus[i]=0;
 109   }
 110
 111   for(Int_t i=0;i<kMaxPtBins+1;i++){
 112     fArrayBinLimits[i]=0;
 113   }
 114
 115 }
 116
 117 //________________________________________________________________________
 118 AliAnalysisTaskSEDplusCorrelations::AliAnalysisTaskSEDplusCorrelations(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana, AliHFAssociatedTrackCuts *asstrkcuts):
 119   AliAnalysisTaskSE(name),
 120   fOutput(0),
 121   fCorrelator(0),
 122   fSelect(0),
 123   fDisplacement(0),
 124   fHistNEvents(0),
 125   fMCSources(0),
 126   fK0Origin(0),
 127   fKaonOrigin(0),
 128   fInvMassK0S(0),
 129   fEventMix(0),
 130   fPtVsMass(0),
 131   fPtVsMassTC(0),
 132   fYVsPt(0),
 133   fYVsPtTC(0),
 134   fYVsPtSig(0),
 135   fYVsPtSigTC(0),
 136   fUpmasslimit(1.965),
 137   fLowmasslimit(1.765),
 138   fNPtBins(0),
 139   fBinWidth(0.002),
 140   fListCuts(0),
 141   fListCutsAsso(0),
 142   fRDCutsAnalysis(dpluscutsana),
 143   fCuts(asstrkcuts),
 144   fCounter(0),
 145   fReadMC(kFALSE),
 146   fUseBit(kTRUE),
 147   fMixing(kFALSE),
 148   fSystem(kFALSE)
 149 {
 150   //
 151   // Standrd constructor
 152   //
 153   fNPtBins=fRDCutsAnalysis->GetNPtBins();
 154
 155   for(Int_t i=0;i<3*kMaxPtBins;i++){
 156     fMassVsdPhiHistHad[i]=0;
 157     fMassVsdEtaHistHad[i]=0;
 158     fMassVsdPhiHistKaon[i]=0;
 159     fMassVsdEtaHistKaon[i]=0;
 160     fMassVsdPhiHistKshort[i]=0;
 161     fMassVsdEtaHistKshort[i]=0;
 162     fMassVsdPhiHistLeadHad[i]=0;
 163     fMassVsdEtaHistLeadHad[i]=0;
 164     fMassHistK0S[i]=0;
 165     fLeadPt[i]=0;
 166     fMassHist[i]=0;
 167     fMassHistTC[i]=0;
 168     fMassHistTCPlus[i]=0;
 169     fMassHistTCMinus[i]=0;
 170   }
 171
 172
 173   for(Int_t i=0;i<kMaxPtBins+1;i++){
 174     fArrayBinLimits[i]=0;
 175   }
 176
 177
 178   // Default constructor
 179   // Output slot #1 writes into a TList container
 180   DefineOutput(1,TList::Class());  //My private output
 181   // Output slot #2 writes cut to private output
 182   //  DefineOutput(2,AliRDHFCutsDplusCorrelationstoKpipi::Class());
 183   DefineOutput(2,TList::Class());
 184   // Output slot #3 writes cut to private output
 185   DefineOutput(3,AliNormalizationCounter::Class());
 186   DefineOutput(4,AliHFAssociatedTrackCuts::Class());
 187
 188
 189 }
 190
 191 //________________________________________________________________________
 192 AliAnalysisTaskSEDplusCorrelations::~AliAnalysisTaskSEDplusCorrelations()
 193 {
 194   //
 195   // Destructor
 196   //
 197   delete fOutput;
 198   delete fListCuts;
 199   delete fRDCutsAnalysis;
 200   delete fCuts;
 201   delete fCounter;
 202   delete fCorrelator;
 203
 204
 205 }
 206 //_________________________________________________________________
 207 void  AliAnalysisTaskSEDplusCorrelations::SetMassLimits(Float_t range){
 208   // set invariant mass limits
 209   Float_t bw=GetBinWidth();
 210   fUpmasslimit = 1.865+range;
 211   fLowmasslimit = 1.865-range;
 212   SetBinWidth(bw);
 213 }
 214 //_________________________________________________________________
 215 void  AliAnalysisTaskSEDplusCorrelations::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 216   // set invariant mass limits
 217   if(uplimit>lowlimit)
 218     {
 219       Float_t bw=GetBinWidth();
 220       fUpmasslimit = lowlimit;
 221       fLowmasslimit = uplimit;
 222       SetBinWidth(bw);
 223     }
 224 }
 225 //________________________________________________________________
 226 void AliAnalysisTaskSEDplusCorrelations::SetBinWidth(Float_t w){
 227   // Define width of mass bins
 228   Float_t width=w;
 229   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
 230   Int_t missingbins=4-nbins%4;
 231   nbins=nbins+missingbins;
 232   width=(fUpmasslimit-fLowmasslimit)/nbins;
 233   if(missingbins!=0){
 234     printf("AliAnalysisTaskSEDplusCorrelations::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
 235   }
 236   else{
 237     if(fDebug>1) printf("AliAnalysisTaskSEDplusCorrelations::SetBinWidth: width set to %f\n",width);
 238   }
 239   fBinWidth=width;
 240 }
 241 //_________________________________________________________________
 242 Int_t AliAnalysisTaskSEDplusCorrelations::GetNBinsHistos(){
 243   // Compute number of mass bins
 244   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
 245 }
 246
 247
 248 //__________________________________________
 249 void AliAnalysisTaskSEDplusCorrelations::Init(){
 250   //
 251   // Initialization
 252   //
 253   if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations::Init() \n");
 254
 255   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
 256   fListCuts=new TList();
 257   // fListCutsAsso=new TList();
 258
 259   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
 260   analysis->SetName("AnalysisCuts");
 261
 262   // AliHFAssociatedTrackCuts *trkcuts = new AliHFAssociatedTrackCuts(*fCuts);
 263   //trkcuts->SetName("Assotrkcuts");
 264
 265   fListCuts->Add(analysis);
 266   //fListCuts->Add(trkcuts);
 267
 268
 269
 270   PostData(2,fListCuts);
 271   PostData(4,fCuts);
 272
 273   return;
 274 }
 275
 276 //________________________________________________________________________
 277 void AliAnalysisTaskSEDplusCorrelations::UserCreateOutputObjects()
 278 {
 279   // Create the output container
 280   //
 281   if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations::UserCreateOutputObjects() \n");
 282   // correlator creation and definition
 283
 284   Double_t Pi = TMath::Pi();
 285   fCorrelator = new AliHFCorrelator("Correlator",fCuts,fSystem); // fCuts is the hadron cut object, fSystem to switch between pp or PbPb
 286   fCorrelator->SetDeltaPhiInterval((-0.5-1./32)*Pi,(1.5-1./32)*Pi); // set correct phi interval
 287   //fCorrelator->SetDeltaPhiInterval(-1.57,4.71);
 288   fCorrelator->SetEventMixing(fMixing); //set kFALSE/kTRUE for mixing Off/On
 289   fCorrelator->SetAssociatedParticleType(fSelect); // set 1/2/3 for hadron/kaons/kzeros
 290   fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
 291   fCorrelator->SetUseMC(fReadMC);
 292   fCorrelator->SetPIDmode(2);
 293
 294   Bool_t pooldef = fCorrelator->DefineEventPool();
 295
 296   if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
 297
 298
 299   // Several histograms are more conveniently managed in a TList
 300   fOutput = new TList();
 301   fOutput->SetOwner();
 302   fOutput->SetName("OutputHistos");
 303
 304   TString hisname;
 305   Int_t index=0;
 306   Int_t nbins=GetNBinsHistos();
 307
 308
 309    Int_t nbinsphi = 32;
 310    Double_t philow = -0.5*Pi - Pi/32; // shift the bin by half the width so that at 0 is it the bin center
 311    Double_t phiup = 1.5*Pi - Pi/32;
 312
 313    Int_t nbinseta = 16;
 314    Double_t etalow = -1.6; // shift the bin by half the width so that at 0 is it the bin center
 315    Double_t etaup = +1.6;
 316
 317
 318
 319   for(Int_t i=0;i<fNPtBins;i++){
 320
 321     index=GetHistoIndex(i);
 322
 323
 324     hisname.Form("hMassVsdPhiHad%d",i);
 325     fMassVsdPhiHistHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 326     fMassVsdPhiHistHad[index]->Sumw2();
 327
 328     hisname.Form("hMassVsdEtaHad%d",i);
 329     fMassVsdEtaHistHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 330     fMassVsdEtaHistHad[index]->Sumw2();
 331
 332     hisname.Form("hMassVsdPhiKaon%d",i);
 333     fMassVsdPhiHistKaon[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 334     fMassVsdPhiHistKaon[index]->Sumw2();
 335
 336     hisname.Form("hMassVsdEtaKaon%d",i);
 337     fMassVsdEtaHistKaon[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 338     fMassVsdEtaHistKaon[index]->Sumw2();
 339
 340     hisname.Form("hMassVsdPhiK0%d",i);
 341     fMassVsdPhiHistKshort[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 342     fMassVsdPhiHistKshort[index]->Sumw2();
 343
 344     hisname.Form("hMassVsdEtaK0%d",i);
 345     fMassVsdEtaHistKshort[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 346     fMassVsdEtaHistKshort[index]->Sumw2();
 347
 348     hisname.Form("hMassVsdPhiLeadHad%d",i);
 349     fMassVsdPhiHistLeadHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 350     fMassVsdPhiHistLeadHad[index]->Sumw2();
 351
 352     hisname.Form("hMassVsdEtaLeadHad%d",i);
 353     fMassVsdEtaHistLeadHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 354     fMassVsdEtaHistLeadHad[index]->Sumw2();
 355
 356     hisname.Form("hMassPtK0S%d",i);
 357     fMassHistK0S[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.3,0.8);
 358     fMassHistK0S[index]->Sumw2();
 359
 360     hisname.Form("hLeadPt%d",i);
 361     fLeadPt[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.0,50);
 362     fLeadPt[index]->Sumw2();
 363
 364
 365     hisname.Form("hMassPt%d",i);
 366     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 367     fMassHist[index]->Sumw2();
 368
 369
 370     hisname.Form("hMassPt%dTC",i);
 371     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 372     fMassHistTC[index]->Sumw2();
 373
 374     hisname.Form("hMassPt%dTCPlus",i);
 375     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 376     fMassHistTCPlus[index]->Sumw2();
 377
 378     hisname.Form("hMassPt%dTCMinus",i);
 379     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 380     fMassHistTCMinus[index]->Sumw2();
 381
 382
 383
 384     index=GetSignalHistoIndex(i);
 385
 386     hisname.Form("hMassVsdPhiHadSig%d",i);
 387     fMassVsdPhiHistHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 388     fMassVsdPhiHistHad[index]->Sumw2();
 389
 390     hisname.Form("hMassVsdEtaHadSig%d",i);
 391     fMassVsdEtaHistHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 392     fMassVsdEtaHistHad[index]->Sumw2();
 393
 394     hisname.Form("hMassVsdPhiKaonSig%d",i);
 395     fMassVsdPhiHistKaon[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 396     fMassVsdPhiHistKaon[index]->Sumw2();
 397
 398     hisname.Form("hMassVsdEtaKaonSig%d",i);
 399     fMassVsdEtaHistKaon[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 400     fMassVsdEtaHistKaon[index]->Sumw2();
 401
 402     hisname.Form("hMassVsdPhiK0Sig%d",i);
 403     fMassVsdPhiHistKshort[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 404     fMassVsdPhiHistKshort[index]->Sumw2();
 405
 406     hisname.Form("hMassVsdEtaK0Sig%d",i);
 407     fMassVsdEtaHistKshort[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 408     fMassVsdEtaHistKshort[index]->Sumw2();
 409
 410     hisname.Form("hMassVsdPhiLeadHadSig%d",i);
 411     fMassVsdPhiHistLeadHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 412     fMassVsdPhiHistLeadHad[index]->Sumw2();
 413
 414     hisname.Form("hMassVsdEtaLeadHadSig%d",i);
 415     fMassVsdEtaHistLeadHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 416     fMassVsdEtaHistLeadHad[index]->Sumw2();
 417
 418
 419     hisname.Form("hSigPtK0S%d",i);
 420     fMassHistK0S[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.3,0.8);
 421     fMassHistK0S[index]->Sumw2();
 422
 423     hisname.Form("hSigLeadPt%d",i);
 424     fLeadPt[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.0,50);
 425     fLeadPt[index]->Sumw2();
 426
 427     hisname.Form("hSigPt%d",i);
 428     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 429     fMassHist[index]->Sumw2();
 430
 431     hisname.Form("hSigPt%dTC",i);
 432     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 433     fMassHistTC[index]->Sumw2();
 434     hisname.Form("hSigPt%dTCPlus",i);
 435     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 436     fMassHistTCPlus[index]->Sumw2();
 437     hisname.Form("hSigPt%dTCMinus",i);
 438     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 439     fMassHistTCMinus[index]->Sumw2();
 440
 441
 442
 443     index=GetBackgroundHistoIndex(i);
 444
 445     hisname.Form("hMassVsdPhiBkgHad%d",i);
 446     fMassVsdPhiHistHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 447     fMassVsdPhiHistHad[index]->Sumw2();
 448
 449     hisname.Form("hMassVsdEtaBkgHad%d",i);
 450     fMassVsdEtaHistHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 451     fMassVsdEtaHistHad[index]->Sumw2();
 452
 453     hisname.Form("hMassVsdPhiBkgKaon%d",i);
 454     fMassVsdPhiHistKaon[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 455     fMassVsdPhiHistKaon[index]->Sumw2();
 456
 457     hisname.Form("hMassVsdEtaBkgKaon%d",i);
 458     fMassVsdEtaHistKaon[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 459     fMassVsdEtaHistKaon[index]->Sumw2();
 460
 461     hisname.Form("hMassVsdPhiBkgKshort%d",i);
 462     fMassVsdPhiHistKshort[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 463     fMassVsdPhiHistKshort[index]->Sumw2();
 464
 465
 466     hisname.Form("hMassVsdPhiBkgKshort%d",i);
 467     fMassVsdEtaHistKshort[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 468     fMassVsdEtaHistKshort[index]->Sumw2();
 469
 470     hisname.Form("hMassVsdPhiBkgLeadHad%d",i);
 471     fMassVsdPhiHistLeadHad[index]=new TH2F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup);
 472     fMassVsdPhiHistLeadHad[index]->Sumw2();
 473
 474     hisname.Form("hMassVsdPhiBkgKshort%d",i);
 475     fMassVsdEtaHistLeadHad[index]=new TH3D(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit,nbinsphi,philow,phiup,nbinseta,etalow,etaup);
 476     fMassVsdEtaHistLeadHad[index]->Sumw2();
 477
 478     hisname.Form("hBkgPtK0S%d",i);
 479     fMassHistK0S[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.3,0.8);
 480     fMassHistK0S[index]->Sumw2();
 481
 482     hisname.Form("hLeadBkgPt%d",i);
 483     fLeadPt[index]=new TH1F(hisname.Data(),hisname.Data(),500,0.0,50);
 484     fLeadPt[index]->Sumw2();
 485
 486     hisname.Form("hBkgPt%dTC",i);
 487     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 488     fMassHistTC[index]->Sumw2();
 489
 490     hisname.Form("hBkgPt%dTCPlus",i);
 491     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 492     fMassHistTCPlus[index]->Sumw2();
 493
 494     hisname.Form("hBkgPt%dTCMinus",i);
 495     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 496     fMassHistTCMinus[index]->Sumw2();
 497   }
 498
 499
 500   for(Int_t i=0; i<3*fNPtBins; i++){
 501     fOutput->Add(fMassVsdPhiHistHad[i]);
 502     fOutput->Add(fMassVsdEtaHistHad[i]);
 503     fOutput->Add(fMassVsdPhiHistKaon[i]);
 504     fOutput->Add(fMassVsdEtaHistKaon[i]);
 505     fOutput->Add(fMassVsdPhiHistKshort[i]);
 506     fOutput->Add(fMassVsdEtaHistKshort[i]);
 507     fOutput->Add(fMassVsdPhiHistLeadHad[i]);
 508     fOutput->Add(fMassVsdEtaHistLeadHad[i]);
 509     fOutput->Add(fMassHistK0S[i]);
 510     fOutput->Add(fLeadPt[i]);
 511     fOutput->Add(fMassHist[i]);
 512     fOutput->Add(fMassHistTC[i]);
 513     fOutput->Add(fMassHistTCPlus[i]);
 514     fOutput->Add(fMassHistTCMinus[i]);
 515
 516
 517
 518   }
 519   fInvMassK0S = new TH2F("K0S","K0S", 500,0.3,0.8,500,0,50);
 520   fInvMassK0S->GetXaxis()->SetTitle("Invariant Mass (#pi #pi) (GeV/c^{2})");
 521   fInvMassK0S->GetYaxis()->SetTitle("K0S pt (GeV/c)");
 522   fOutput->Add(fInvMassK0S);
 523
 524   fMCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
 525   fMCSources->GetXaxis()->SetBinLabel(1,"All ");
 526   fMCSources->GetXaxis()->SetBinLabel(2," from Heavy flavour");
 527   fMCSources->GetXaxis()->SetBinLabel(3," from c->D");
 528   fMCSources->GetXaxis()->SetBinLabel(4," from b->D");
 529   fMCSources->GetXaxis()->SetBinLabel(5," from b->B");
 530   if(fReadMC) fOutput->Add(fMCSources);
 531
 532   fK0Origin = new TH1F("K0Origin","Origin of K0 ", 10, -0.5, 5.5);
 533   fK0Origin->GetXaxis()->SetBinLabel(1,"All K0s");
 534   fK0Origin->GetXaxis()->SetBinLabel(2,"K0s from Heavy flavour");
 535   fK0Origin->GetXaxis()->SetBinLabel(3,"K0s from Charm");
 536   fK0Origin->GetXaxis()->SetBinLabel(4,"K0s from Beauty");
 537   if(fReadMC) fOutput->Add(fK0Origin);
 538
 539   fKaonOrigin = new TH1F("K0Origin","Origin of Kaon ", 10, -0.5, 5.5);
 540   fKaonOrigin->GetXaxis()->SetBinLabel(1,"All Kaons");
 541   fKaonOrigin->GetXaxis()->SetBinLabel(2,"Kaons from Heavy flavour");
 542   fKaonOrigin->GetXaxis()->SetBinLabel(3,"Kaons from Charm");
 543   fKaonOrigin->GetXaxis()->SetBinLabel(4,"Kaons from Beauty");
 544   if(fReadMC) fOutput->Add(fKaonOrigin);
 545
 546
 547
 548   fEventMix = new TH2F("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
 549   if(fMixing)fOutput->Add(fEventMix);
 550
 551
 552   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 553   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 554   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
 555   fHistNEvents->GetXaxis()->SetBinLabel(3,"Rejected due to trigger");
 556   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected pileup events");
 557   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to centrality");
 558   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vtxz");
 559   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to Physics Sel");
 560   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 561   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 562   fHistNEvents->GetXaxis()->SetBinLabel(10,"D+ after loose cuts");
 563   fHistNEvents->GetXaxis()->SetBinLabel(11,"D+ after tight cuts");
 564
 565   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 566   fHistNEvents->Sumw2();
 567   fHistNEvents->SetMinimum(0);
 568   fOutput->Add(fHistNEvents);
 569
 570   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 571   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 572   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 573   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
 574   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 575   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 576
 577   fOutput->Add(fPtVsMass);
 578   fOutput->Add(fPtVsMassTC);
 579   fOutput->Add(fYVsPt);
 580   fOutput->Add(fYVsPtTC);
 581   fOutput->Add(fYVsPtSig);
 582   fOutput->Add(fYVsPtSigTC);
 583
 584
 585   // Counter for Normalization
 586   TString normName="NormalizationCounter";
 587   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
 588   if(cont)normName=(TString)cont->GetName();
 589   fCounter = new AliNormalizationCounter(normName.Data());
 590   fCounter->Init();
 591
 592
 593
 594   PostData(1,fOutput);
 595   PostData(3,fCounter);
 596   return;
 597 }
 598
 599 //________________________________________________________________________
 600 void AliAnalysisTaskSEDplusCorrelations::UserExec(Option_t */*option*/)
 601 {
 602   // Do the analysis
 603   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 604   TClonesArray *array3Prong = 0;
 605
 606   if(!fMixing){
 607     if(fSelect==1) cout << "TASK::Correlation with hadrons on SE "<< endl;
 608     if(fSelect==2) cout << "TASK::Correlation with kaons on SE "<< endl;
 609     if(fSelect==3) cout << "TASK::Correlation with kzeros on SE "<< endl;
 610   }
 611   if(fMixing){
 612     if(fSelect==1) cout << "TASK::Correlation with hadrons on ME "<< endl;
 613     if(fSelect==2) cout << "TASK::Correlation with kaons on ME "<< endl;
 614     if(fSelect==3) cout << "TASK::Correlation with kzeros on ME "<< endl;
 615   }
 616
 617
 618   if(!aod && AODEvent() && IsStandardAOD()) {
 619
 620     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 621     AliAODHandler* aodHandler = (AliAODHandler*)
 622       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 623     if(aodHandler->GetExtensions()) {
 624       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 625       AliAODEvent *aodFromExt = ext->GetAOD();
 626       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 627     }
 628   } else if(aod) {
 629     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 630   }
 631
 632   if(!aod || !array3Prong) {
 633     printf("AliAnalysisTaskSEDplusCorrelations::UserExec: Charm3Prong branch not found!\n");
 634     return;
 635   }
 636
 637
 638   // the AODs with null vertex pointer didn't pass the PhysSel
 639   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 640   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
 641   fHistNEvents->Fill(0); // count event
 642
 643
 644   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 645
 646
 647
 648   PostData(1,fOutput);
 649   if(!isEvSel)return;
 650   fHistNEvents->Fill(1);
 651
 652   // set PIDResponse for associated tracks
 653   fCorrelator->SetPidAssociated();
 654
 655   TClonesArray *arrayMC=0;
 656   AliAODMCHeader *mcHeader=0;
 657   // AOD primary vertex
 658   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 659   //    vtx1->Print();
 660   TString primTitle = vtx1->GetTitle();
 661   //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 662
 663   // load MC particles
 664   if(fReadMC){
 665
 666     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 667     if(!arrayMC) {
 668       printf("AliAnalysisTaskSEDplusCorrelations::UserExec: MC particles branch not found!\n");
 669       return;
 670     }
 671
 672     // load MC header
 673     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 674     if(!mcHeader) {
 675       printf("AliAnalysisTaskSEDplusCorrelations::UserExec: MC header branch not found!\n");
 676       return;
 677     }
 678   }
 679
 680   //HFCorrelators initialization (for this event)
 681
 682   fCorrelator->SetAODEvent(aod); // set the event to be processedfCorrelator->
 683   Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
 684   if(!correlatorON) {
 685     AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
 686     return;
 687   }
 688   if(fReadMC) fCorrelator->SetMCArray(arrayMC);
 689
 690
 691
 692   Int_t n3Prong = array3Prong->GetEntriesFast();
 693    printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
 694   Int_t nOS=0;
 695   Int_t index;
 696   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 697   Int_t nSelectedloose=0,nSelectedtight=0;
 698
 699
 700
 701   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 702     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 703     fHistNEvents->Fill(7);
 704
 705     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 706       fHistNEvents->Fill(8);
 707       continue;
 708     }
 709
 710     Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 711
 712      if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
 713
 714     Bool_t unsetvtx=kFALSE;
 715     if(!d->GetOwnPrimaryVtx()){
 716       d->SetOwnPrimaryVtx(vtx1);
 717       unsetvtx=kTRUE;
 718     }
 719
 720
 721     Double_t ptCand = d->Pt();
 722     Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
 723     Bool_t recVtx=kFALSE;
 724     AliAODVertex *origownvtx=0x0;
 725     if(fRDCutsAnalysis->GetIsPrimaryWithoutDaughters()){
 726       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 727       if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
 728       else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
 729     }
 730
 731
 732
 733
 734     Int_t labDp=-1;
 735     Bool_t isDplus = kFALSE;
 736
 737     if(fReadMC){
 738       labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
 739       if(labDp>=0){
 740         isDplus = kTRUE;
 741       }
 742     }
 743
 744
 745     Double_t invMass=d->InvMassDplus();
 746     Double_t rapid=d->YDplus();
 747     fYVsPt->Fill(ptCand,rapid);
 748     if(passTightCuts>0.) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
 749     //printf("****************: %d and of tracks: %d\n",n3Prong,nOS);
 750     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 751     if(isFidAcc){
 752       fPtVsMass->Fill(invMass,ptCand);
 753       if(passTightCuts>0) fPtVsMassTC->Fill(invMass,ptCand);
 754     }
 755
 756
 757     if(iPtBin>=0){
 758       index=GetHistoIndex(iPtBin);
 759       if(isFidAcc){
 760         fHistNEvents->Fill(9);
 761         nSelectedloose++;
 762         fMassHist[index]->Fill(invMass);
 763
 764         // loop for tight cuts
 765         if(passTightCuts>0){
 766           fHistNEvents->Fill(10);
 767           nSelectedtight++;
 768           fMassHistTC[index]->Fill(invMass);
 769
 770           //Dplus info
 771
 772           Double_t phiDplus = fCorrelator->SetCorrectPhiRange(d->Phi());
 773           fCorrelator->SetTriggerParticleProperties(d->Pt(),phiDplus,d->Eta());
 774
 775           Int_t nIDs[3] = {-9999999};
 776           nIDs[0] = ((AliAODTrack*)d->GetDaughter(0))->GetID();
 777           nIDs[1] = ((AliAODTrack*)d->GetDaughter(1))->GetID();
 778           nIDs[2] = ((AliAODTrack*)d->GetDaughter(2))->GetID();
 779
 780           Double_t ptlead = 0;
 781           Double_t philead = 0;
 782           Double_t etalead = 0;
 783           Double_t refpt = 0;
 784
 785
 786
 787           Bool_t execPool = fCorrelator->ProcessEventPool();
 788
 789           //     printf("*************: %d\n",execPool);
 790           if(fMixing && !execPool) {
 791             AliInfo("Mixed event analysis: pool is not ready");
 792             continue;
 793           }
 794           Int_t NofEventsinPool = 1;
 795           if(fMixing) {
 796             NofEventsinPool = fCorrelator->GetNofEventsInPool();
 797           }
 798
 799
 800
 801           for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
 802
 803             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
 804             if(!analyzetracks) {
 805               AliInfo("AliHFCorrelator::Cannot process the track array");
 806               continue;
 807             }
 808
 809             //start the track loop
 810
 811             // Int_t NofTracks = fCorrelator->GetNofTracks();
 812
 813             //cout<<"*******"<<NofTracks<<endl;
 814
 815             for (Int_t iTrack = 0;iTrack<fCorrelator->GetNofTracks();iTrack++) {
 816               Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
 817
 818               if(!runcorrelation) continue;
 819               Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
 820               Double_t DeltaEta = fCorrelator->GetDeltaEta();
 821
 822               AliReducedParticle* redpart = fCorrelator->GetAssociatedParticle();
 823               Double_t phiHad=redpart->Phi();
 824               Double_t ptHad=redpart->Pt();
 825               Double_t etaHad=redpart->Eta();
 826               Int_t label = redpart->GetLabel();
 827               Int_t trackid = redpart->GetID();
 828               phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
 829
 830
 831               //  discard the dplus daughters
 832               if (!fMixing){
 833                 if( trackid == nIDs[0] || trackid == nIDs[1] || trackid == nIDs[2]) continue;
 834               }
 835               // discard the negative id tracks
 836               if(trackid < 0) continue;
 837
 838
 839               FillCorrelations(d,DeltaPhi,DeltaEta,index,fSelect);
 840
 841               // For leading particle
 842
 843               if (ptHad > refpt) {
 844                 refpt = ptHad; ptlead = ptHad;
 845                 philead = d->Phi() - phiHad;
 846                 etalead = d->Eta() - etaHad;
 847                 if (philead < (-1)*TMath::Pi()/2) philead += 2*TMath::Pi();
 848                 if (philead > 3*TMath::Pi()/2) philead -= 2*TMath::Pi();
 849
 850               }
 851
 852               // montecarlo
 853
 854               if(fReadMC && isDplus) {
 855
 856                 index=GetSignalHistoIndex(iPtBin);
 857
 858                 Bool_t* partSource = fCuts->IsMCpartFromHF(label,arrayMC); // check source of associated particle (hadron/kaon/K0)
 859                 FillMCCorrelations(d,DeltaPhi,DeltaEta,index,partSource,fSelect);
 860                 delete partSource;
 861
 862
 863               } // readMC
 864
 865             }//count good tracks
 866
 867             // For leading particle
 868             fMassVsdPhiHistLeadHad[index]->Fill(invMass,philead);
 869             fMassVsdEtaHistLeadHad[index]->Fill(invMass,philead,etalead);
 870
 871             fLeadPt[index]->Fill(ptlead);
 872
 873             if(fReadMC && isDplus) {
 874               index=GetSignalHistoIndex(iPtBin);
 875               fMassVsdPhiHistLeadHad[index]->Fill(invMass,philead);
 876               fMassVsdEtaHistLeadHad[index]->Fill(invMass,philead,etalead);
 877               fLeadPt[index]->Fill(ptlead);
 878
 879             }
 880
 881           }//jmix
 882
 883           }// tc
 884       }//fid
 885
 886
 887
 888
 889
 890     }
 891     if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
 892
 893     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 894
 895   }
 896
 897   if(fMixing){
 898     Bool_t updated = fCorrelator->PoolUpdate();
 899
 900     if(!updated) AliInfo("Pool was not updated");
 901   }
 902   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 903   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 904
 905
 906   PostData(1,fOutput);
 907   PostData(2,fListCuts);
 908   PostData(3,fCounter);
 909   return;
 910 }
 911
 912 //________________________________________________________________________
 913 void AliAnalysisTaskSEDplusCorrelations::Terminate(Option_t */*option*/)
 914 {
 915   // Terminate analysis
 916   //
 917   if(fDebug > 1) printf("AnalysisTaskSEDplusCorrelations: Terminate() \n");
 918
 919   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 920   if (!fOutput) {
 921     printf("ERROR: fOutput not available\n");
 922     return;
 923   }
 924   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 925
 926   TString hisname;
 927   Int_t index=0;
 928
 929   for(Int_t i=0;i<fNPtBins;i++){
 930     index=GetHistoIndex(i);
 931
 932     hisname.Form("hMassPt%dTC",i);
 933     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 934   }
 935
 936   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
 937   c1->cd();
 938   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
 939   hMassPt->SetLineColor(kBlue);
 940   hMassPt->SetXTitle("M[GeV/c^{2}]");
 941   hMassPt->Draw();
 942
 943   return;
 944 }
 945
 946
 947 //________________________________________________________________________
 948 void AliAnalysisTaskSEDplusCorrelations::FillCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t sel) const{
 949   //Filling histogams
 950
 951   Double_t invMass=d->InvMassDplus();
 952
 953
 954   if(sel==1){
 955     fMassVsdPhiHistHad[ind]->Fill(invMass,deltaPhi);
 956     fMassVsdEtaHistHad[ind]->Fill(invMass,deltaPhi,deltaEta);
 957   }
 958   if(sel==2){
 959     fMassVsdPhiHistKaon[ind]->Fill(invMass,deltaPhi);
 960     fMassVsdEtaHistKaon[ind]->Fill(invMass,deltaPhi,deltaEta);
 961   }
 962   if(sel==3){
 963     fMassVsdPhiHistKshort[ind]->Fill(invMass,deltaPhi);
 964     fMassVsdEtaHistKshort[ind]->Fill(invMass,deltaPhi,deltaEta);
 965   }
 966
 967   return;
 968 }
 969
 970
 971
 972 //________________________________________________________________________
 973 void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind,Bool_t* mcSource, Int_t sel) const{
 974   // Filling histos with MC information
 975
 976   Double_t invMass=d->InvMassDplus();
 977
 978
 979   if(sel==1){
 980     fMassVsdPhiHistHad[ind]->Fill(invMass,deltaPhi);
 981     fMassVsdEtaHistHad[ind]->Fill(invMass,deltaPhi,deltaEta);
 982
 983
 984     fMCSources->Fill(0);
 985
 986     if(mcSource[2]&&mcSource[0]){ // is from charm ->D
 987       fMCSources->Fill(1);
 988       fMCSources->Fill(2);
 989     }
 990         if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
 991       fMCSources->Fill(1);
 992       fMCSources->Fill(3);
 993     }
 994     if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
 995       fMCSources->Fill(1);
 996       fMCSources->Fill(4);
 997     }
 998   }
 999
1000   if(sel==2){
1001     fMassVsdPhiHistKaon[ind]->Fill(invMass,deltaPhi);
1002     fMassVsdEtaHistKaon[ind]->Fill(invMass,deltaPhi,deltaEta);
1003     fKaonOrigin->Fill(0);
1004     if(mcSource[2]&&mcSource[0]){ // is from charm ->D
1005       fKaonOrigin->Fill(1);
1006       fKaonOrigin->Fill(2);
1007     }
1008     if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
1009       fKaonOrigin->Fill(1);
1010       fKaonOrigin->Fill(3);
1011     }
1012     if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
1013       fKaonOrigin->Fill(1);
1014       fKaonOrigin->Fill(4);
1015     }
1016   }
1017   if(sel==3){
1018     fMassVsdPhiHistKshort[ind]->Fill(invMass,deltaPhi);
1019     fMassVsdEtaHistKshort[ind]->Fill(invMass,deltaPhi,deltaEta);
1020     fK0Origin->Fill(0);
1021     if(mcSource[2]&&mcSource[0]){ // is from charm ->D
1022       fK0Origin->Fill(1);
1023       fK0Origin->Fill(2);
1024     }
1025     if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
1026       fK0Origin->Fill(1);
1027       fK0Origin->Fill(3);
1028     }
1029     if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
1030       fK0Origin->Fill(1);
1031       fK0Origin->Fill(4);
1032     }
1033
1034   }
1035
1036   return;
1037 }
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047