PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 //
  17 // Task for selecting D mesons to be used as an input for D-Jet correlations
  18 //
  19 //-----------------------------------------------------------------------
  20 // Authors:
  21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
  22 // A.Grelli (Utrecht University) a.grelli@uu.nl
  23 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
  24 //-----------------------------------------------------------------------
  25
  26 #include <TDatabasePDG.h>
  27 #include <TParticle.h>
  28 #include <TVector3.h>
  29 #include "TROOT.h"
  30 #include <TH3F.h>
  31
  32 #include "AliRDHFCutsDStartoKpipi.h"
  33 #include "AliRDHFCutsD0toKpi.h"
  34 #include "AliAODMCHeader.h"
  35 #include "AliAODHandler.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliLog.h"
  38 #include "AliAODVertex.h"
  39 #include "AliAODRecoDecay.h"
  40 #include "AliAODRecoCascadeHF.h"
  41 #include "AliAODRecoDecayHF2Prong.h"
  42 #include "AliESDtrack.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
  45
  46 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
  47
  48 //_______________________________________________________________________________
  49
  50 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
  51 AliAnalysisTaskSE(),
  52 fUseMCInfo(kFALSE),
  53 fUseReco(kTRUE),
  54 fCandidateType(0),
  55 fCandidateName(""),
  56 fPDGmother(0),
  57 fNProngs(0),
  58 fBranchName(""),
  59 fOutput(0),
  60 fCuts(0),
  61 fMinMass(0.),
  62 fMaxMass(0.),
  63 fCandidateArray(0),
  64 fSideBandArray(0),
  65 fhImpPar(0),
  66 fhImpParB(0),
  67 fhInvMassS(0),
  68 fhInvMassB(0)
  69 {
  70    //
  71    // Default constructor
  72    //
  73
  74    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
  75    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
  76 }
  77
  78 //_______________________________________________________________________________
  79
  80 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
  81 AliAnalysisTaskSE(name),
  82 fUseMCInfo(kFALSE),
  83 fUseReco(kTRUE),
  84 fCandidateType(candtype),
  85 fCandidateName(""),
  86 fPDGmother(0),
  87 fNProngs(0),
  88 fBranchName(""),
  89 fOutput(0),
  90 fCuts(cuts),
  91 fMinMass(0.),
  92 fMaxMass(0.),
  93 fCandidateArray(0),
  94 fSideBandArray(0),
  95 fhImpPar(0),
  96 fhImpParB(0),
  97 fhInvMassS(0),
  98 fhInvMassB(0)
  99 {
 100    //
 101    // Constructor. Initialization of Inputs and Outputs
 102    //
 103
 104    Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
 105
 106    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
 107    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
 108
 109    const Int_t nptbins = fCuts->GetNPtBins();
 110    Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
 111
 112    switch (fCandidateType) {
 113    case 0 :
 114       fCandidateName = "D0";
 115       fPDGmother = 421;
 116       fNProngs = 2;
 117       fPDGdaughters[0] = 211;  // pi
 118       fPDGdaughters[1] = 321;  // K
 119       fPDGdaughters[2] = 0;    // empty
 120       fPDGdaughters[3] = 0;    // empty
 121       fBranchName = "D0toKpi";
 122       break;
 123    case 1 :
 124       fCandidateName = "Dstar";
 125       fPDGmother = 413;
 126       fNProngs = 3;
 127       fPDGdaughters[1] = 211; // pi soft
 128       fPDGdaughters[0] = 421; // D0
 129       fPDGdaughters[2] = 211; // pi fromD0
 130       fPDGdaughters[3] = 321; // K from D0
 131       fBranchName = "Dstar";
 132
 133       if (nptbins<=13) {
 134          for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
 135       } else {
 136          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 137       }
 138       break;
 139    default :
 140       printf("%d not accepted!!\n",fCandidateType);
 141       break;
 142    }
 143
 144    if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
 145    if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 146
 147    DefineOutput(1, TList::Class());       // histos
 148    DefineOutput(2, AliRDHFCuts::Class()); // my cuts
 149    DefineOutput(3, TClonesArray::Class()); //array of candidates
 150    DefineOutput(4, TClonesArray::Class()); //array of SB candidates
 151 }
 152
 153 //_______________________________________________________________________________
 154
 155 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
 156 {
 157    //
 158    // destructor
 159    //
 160
 161    Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
 162
 163    if (fOutput) { delete fOutput; fOutput = 0; }
 164    if (fCuts)   { delete fCuts;   fCuts   = 0; }
 165    if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
 166    delete fSideBandArray;
 167
 168 }
 169
 170 //_______________________________________________________________________________
 171
 172 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
 173 {
 174    //
 175    // Initialization
 176    //
 177
 178    if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
 179
 180    switch (fCandidateType) {
 181    case 0:
 182       {
 183          AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 184          copyfCutsDzero->SetName("AnalysisCutsDzero");
 185          PostData(2, copyfCutsDzero);  // Post the data
 186       } break;
 187    case 1:
 188       {
 189          AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 190          copyfCutsDstar->SetName("AnalysisCutsDStar");
 191          PostData(2, copyfCutsDstar); // Post the cuts
 192       } break;
 193    default: return;
 194    }
 195
 196    return;
 197 }
 198
 199 //_______________________________________________________________________________
 200
 201 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
 202 {
 203    //
 204    // output
 205    //
 206
 207    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 208
 209    fOutput = new TList(); fOutput->SetOwner();
 210    DefineHistoForAnalysis(); // define histograms
 211
 212    if (fCandidateType==kD0toKpi){
 213       fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
 214       fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0);
 215    }
 216
 217    if (fCandidateType==kDstartoKpipi) {
 218       fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
 219       fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0);
 220    }
 221
 222    fCandidateArray->SetOwner();
 223    fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
 224
 225    //this is used for the DStar side bands and MC!
 226    fSideBandArray->SetOwner();
 227    fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
 228
 229    PostData(1, fOutput);
 230    PostData(3, fCandidateArray);
 231    PostData(4, fSideBandArray);
 232
 233    return;
 234 }
 235
 236 //_______________________________________________________________________________
 237
 238 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
 239    //
 240    // user exec
 241    //
 242    // Load the event
 243    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 244
 245    TClonesArray *arrayDStartoD0pi = 0;
 246    if (!aodEvent && AODEvent() && IsStandardAOD()) {
 247
 248       // In case there is an AOD handler writing a standard AOD, use the AOD
 249       // event in memory rather than the input (ESD) event.
 250       aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
 251
 252       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 253       // have to taken from the AOD event hold by the AliAODExtension
 254       AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 255       if(aodHandler->GetExtensions()) {
 256          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 257          AliAODEvent *aodFromExt = ext->GetAOD();
 258          arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 259       }
 260    } else {
 261       arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 262    }
 263
 264    if (!arrayDStartoD0pi) {
 265       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 266       return;
 267    } else {
 268       AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 269    }
 270
 271    TClonesArray* mcArray = 0x0;
 272    if (fUseMCInfo) {
 273       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 274       if (!mcArray) {
 275          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 276          return;
 277       }
 278    }
 279
 280    //Histograms
 281    TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
 282    TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
 283    TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
 284    TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
 285
 286    TH1F* hPtPion=0x0;
 287    if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
 288    hstat->Fill(0);
 289
 290    // fix for temporary bug in ESDfilter
 291    // the AODs with null vertex pointer didn't pass the PhysSel
 292    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 293
 294    //Event selection
 295    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 296    //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 297    if (!iseventselected) return;
 298    hstat->Fill(1);
 299
 300    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
 301    if(!fUseMCInfo) hstat->Fill(2, nD);
 302
 303    //D* and D0 prongs needed to MatchToMC method
 304    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
 305    Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
 306
 307    //D0 from D0 bar
 308    Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K
 309    Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi
 310
 311    Int_t iCand =0;
 312    Int_t iSBCand=0;
 313    Int_t isSelected = 0;
 314    AliAODRecoDecayHF *charmCand = 0;
 315    AliAODRecoCascadeHF *dstar = 0;
 316    AliAODMCParticle *charmPart = 0;
 317    Bool_t isMCBkg=kFALSE;
 318
 319    Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 320
 321    Int_t mcLabel = -9999;
 322    Int_t pdgMeson = 413;
 323    if (fCandidateType==kD0toKpi) pdgMeson = 421;
 324
 325    //clear the TClonesArray from the previous event
 326    fCandidateArray->Clear();
 327    fSideBandArray->Clear();
 328
 329    for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
 330       charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
 331       if (!charmCand) continue;
 332
 333       TString smcTruth="S";
 334
 335       if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
 336
 337       if (fUseMCInfo) { // Look in MC, try to simulate the z
 338          if (fCandidateType==kDstartoKpipi) {
 339             mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 340          }
 341
 342          if (fCandidateType==kD0toKpi)
 343             mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
 344
 345          if (mcLabel<=0) {
 346             isMCBkg=kTRUE;
 347             hstat->Fill(5);
 348          }
 349          else {
 350             hstat->Fill(2);
 351
 352          }
 353          if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
 354
 355          if (isMCBkg) smcTruth="B";
 356
 357       }
 358
 359       Double_t ptD = charmCand->Pt();
 360
 361       // region of interest + cuts
 362       if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
 363
 364       if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
 365          //select by track cuts the side band candidates (don't want mass cut)
 366          isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
 367          if (!isSelected) continue;
 368          //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
 369          Int_t bin = fCuts->PtBin(ptD);
 370          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 371             AliError(Form("Pt %.3f out of bounds",ptD));
 372             continue;
 373          }
 374          //if data and Dstar from D0 side band
 375          if (((dstar->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((dstar->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
 376
 377             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
 378             iSBCand++;
 379          }
 380       }
 381       //candidate selected by cuts and PID
 382       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
 383       if (!isSelected) continue;
 384
 385       Int_t nprongs=charmCand->GetNProngs();
 386
 387       //for data and MC signal fill fCandidateArray
 388       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
 389          // for data or MC with the requirement fUseReco fill with candidates
 390
 391          if(fUseReco) {
 392             if (fCandidateType==kDstartoKpipi){
 393                new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
 394                AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
 395                fhInvMassS->Fill(dstar->DeltaInvMass());
 396
 397             } else{
 398                new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
 399                //Printf("Filling reco");
 400                for(Int_t kd=0;kd<nprongs;kd++){
 401                   fhImpPar->Fill(charmCand->Getd0Prong(kd),charmCand->PtProng(kd));
 402                }
 403                Double_t masses[2];
 404                masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 405                masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 406                fhInvMassS->Fill(masses[0]);
 407                fhInvMassS->Fill(masses[1]);
 408
 409             }
 410
 411             hstat->Fill(3);
 412
 413          }
 414          // for MC with requirement particle level fill with AliAODMCParticle
 415          else if (fUseMCInfo) {
 416             new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
 417             //Printf("Filling gen");
 418             hstat->Fill(3);
 419          }
 420
 421          iCand++;
 422       }
 423       //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
 424       else if(fUseReco){
 425          if (fCandidateType==kDstartoKpipi){
 426             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
 427             fhInvMassB->Fill(dstar->DeltaInvMass());
 428          }
 429          if (fCandidateType==kD0toKpi){
 430             new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
 431             Double_t masses[2];
 432             masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 433             masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 434             fhInvMassB->Fill(masses[0]);
 435             fhInvMassB->Fill(masses[1]);
 436
 437          }
 438          iSBCand++;
 439       }
 440
 441
 442       Double_t masses[2];
 443       if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
 444
 445          //softpion from D* decay
 446          AliAODTrack *track2 = (AliAODTrack*)dstar->GetBachelor();
 447
 448          // select D* in the D0 window.
 449          // In the cut object window is loose to allow for side bands
 450
 451
 452          // retrieve the corresponding pt bin for the candidate
 453          // and set the expected D0 width (x3)
 454          //    static const Int_t n = fCuts->GetNPtBins();
 455          Int_t bin = fCuts->PtBin(ptD);
 456          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 457             AliError(Form("Pt %.3f out of bounds",ptD));
 458             continue;
 459          }
 460
 461          AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
 462          //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
 463          if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
 464             masses[0] = dstar->DeltaInvMass(); //D*
 465             masses[1] = 0.; //dummy for D*
 466
 467             //D*  delta mass
 468             hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
 469
 470             // D* pt and soft pion pt for good candidates
 471             Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 472             Double_t invmassDelta = dstar->DeltaInvMass();
 473             if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
 474          }
 475
 476          if (fUseMCInfo){ //fill histograms of kinematics, using MC truth
 477             //get histos
 478             TH2F *halphaDD   = (TH2F*)fOutput->FindObject(Form("halphaDD%s",smcTruth.Data()));
 479             TH2F *halphaDpis = (TH2F*)fOutput->FindObject(Form("halphaDpis%s",smcTruth.Data()));
 480             TH2F *halphaDpi  = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
 481             TH2F *halphaDK   = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
 482
 483             TH2F *hdeltaRDD   = (TH2F*)fOutput->FindObject(Form("hdeltaRDD%s",smcTruth.Data()));
 484             TH2F *hdeltaRDpis = (TH2F*)fOutput->FindObject(Form("hdeltaRDpis%s",smcTruth.Data()));
 485             TH2F *hdeltaRDpi  = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
 486             TH2F *hdeltaRDK   = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
 487
 488             Double_t aD  = dstar->Phi(),
 489                      apis= track2->Phi();
 490
 491             AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
 492             Double_t aD0 = D0fromDstar->Phi();
 493             Int_t isD0= D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
 494             Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1),
 495                      api= isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
 496             Double_t dRDD0  = DeltaR(dstar,D0fromDstar),
 497                      dRDpis = DeltaR(dstar,track2),
 498                      dRDpi  = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(1) : (AliVParticle*)D0fromDstar->GetDaughter(0)),
 499                      dRDK   = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(0) : (AliVParticle*)D0fromDstar->GetDaughter(1));
 500
 501             halphaDD->  Fill(aD-aD0,ptD);
 502             halphaDpis->Fill(aD-apis,ptD);
 503             halphaDpi-> Fill(aD-api,ptD);
 504             halphaDK->  Fill(aD-aK,ptD);
 505
 506             hdeltaRDD->  Fill(dRDD0,ptD);
 507             hdeltaRDpis->Fill(dRDpis,ptD);
 508             hdeltaRDpi-> Fill(dRDpi,ptD);
 509             hdeltaRDK->  Fill(dRDK,ptD);
 510
 511             if (isMCBkg){
 512                fhImpParB->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
 513                fhImpParB->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
 514                fhImpParB->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
 515
 516             } else {
 517                fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
 518                fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
 519                fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
 520
 521             }
 522
 523
 524          } else{
 525             fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
 526             AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
 527             fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
 528             fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
 529          }
 530
 531       } //Dstar specific
 532
 533       if (fCandidateType==kD0toKpi) { //D0->Kpi
 534
 535          //needed quantities
 536          masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 537          masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 538          hstat->Fill(3);
 539
 540          // mass vs pt
 541          if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
 542          if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
 543
 544          if (fUseMCInfo) {  //fill histograms of kinematics, using MC truth
 545
 546             Double_t aD = charmCand->Phi();
 547             Double_t adaugh[2]={charmCand->PhiProng(0),charmCand->PhiProng(1)};
 548             AliAODTrack* p0=(AliAODTrack*)charmCand->GetDaughter(0);
 549             AliAODTrack* p1=(AliAODTrack*)charmCand->GetDaughter(1);
 550             Float_t dR0 = DeltaR(charmCand, p0), dR1 = DeltaR(charmCand, p1);
 551             Bool_t isD0=kFALSE;
 552             if(mcLabel==421)  isD0=kTRUE;
 553             if(mcLabel==-421) isD0=kFALSE;
 554
 555             if(isMCBkg) { //background
 556                TH2F *halphaDpi  = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
 557                TH2F *halphaDK   = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
 558
 559                TH2F *hdeltaRDpi  = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
 560                TH2F *hdeltaRDK   = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
 561
 562
 563                if (isSelected==1 || isSelected==3) { // selected as D0
 564                   halphaDK->Fill(aD-adaugh[0],ptD);
 565                   halphaDpi->Fill(aD-adaugh[1],ptD);
 566
 567                   hdeltaRDK->Fill(dR0,ptD);
 568                   hdeltaRDpi->Fill(dR1,ptD);
 569
 570                }
 571                if (isSelected>=2) { //selected as D0bar
 572                   halphaDpi->Fill(aD-adaugh[0],ptD);
 573                   halphaDK->Fill(aD-adaugh[1],ptD);
 574
 575                   hdeltaRDpi->Fill(dR0,ptD);
 576                   hdeltaRDK->Fill(dR1,ptD);
 577
 578                }
 579
 580             }else{ //signal and reflections
 581                TH2F *halphaDpiS  = (TH2F*)fOutput->FindObject("halphaDpiS");
 582                TH2F *halphaDKS   = (TH2F*)fOutput->FindObject("halphaDKS");
 583                TH2F *halphaDpiR  = (TH2F*)fOutput->FindObject("halphaDpiR");
 584                TH2F *halphaDKR   = (TH2F*)fOutput->FindObject("halphaDKR");
 585
 586                TH2F *hdeltaRDpiS  = (TH2F*)fOutput->FindObject("hdeltaRDpiS");
 587                TH2F *hdeltaRDKS   = (TH2F*)fOutput->FindObject("hdeltaRDKS");
 588                TH2F *hdeltaRDpiR  = (TH2F*)fOutput->FindObject("hdeltaRDpiR");
 589                TH2F *hdeltaRDKR   = (TH2F*)fOutput->FindObject("hdeltaRDKR");
 590
 591                if(isD0) { //D0
 592                   halphaDKS->Fill(aD-adaugh[0],ptD);
 593                   halphaDpiS->Fill(aD-adaugh[1],ptD);
 594
 595                   hdeltaRDKS->Fill(dR0,ptD);
 596                   hdeltaRDpiS->Fill(dR1,ptD);
 597                   if(isSelected>=2){ //selected as D0bar
 598                      halphaDpiR->Fill(aD-adaugh[0],ptD);
 599                      halphaDKR->Fill(aD-adaugh[1],ptD);
 600
 601                      hdeltaRDpiR->Fill(dR0,ptD);
 602                      hdeltaRDKR->Fill(dR1,ptD);
 603                   }
 604                } else { //D0bar
 605                   halphaDKS->Fill(aD-adaugh[1],ptD);
 606                   halphaDpiS->Fill(aD-adaugh[0],ptD);
 607
 608                   hdeltaRDKS->Fill(dR1,ptD);
 609                   hdeltaRDpiS->Fill(dR0,ptD);
 610
 611                   if(isSelected>=2){ //selected as D0bar
 612                      halphaDpiR->Fill(aD-adaugh[1],ptD);
 613                      halphaDKR->Fill(aD-adaugh[0],ptD);
 614
 615                      hdeltaRDpiR->Fill(dR1,ptD);
 616                      hdeltaRDKR->Fill(dR0,ptD);
 617                   }
 618                }
 619
 620             } //end signal and reflections
 621
 622
 623          }// end MC
 624
 625
 626       } //D0 specific
 627
 628       charmCand = 0;
 629       isMCBkg=kFALSE;
 630    } // end of D cand loop
 631
 632    hnCandEv->Fill(fCandidateArray->GetEntriesFast());
 633    if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
 634       Int_t nsbcand=fSideBandArray->GetEntriesFast();
 635       hstat->Fill(4,nsbcand);
 636       hnSBCandEv->Fill(nsbcand);
 637    }
 638    //Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
 639    PostData(1, fOutput);
 640    PostData(3, fCandidateArray);
 641    PostData(4, fSideBandArray);
 642
 643    return;
 644 }
 645
 646 //_______________________________________________________________________________
 647
 648 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
 649 {
 650    // The Terminate() function is the last function to be called during
 651    // a query. It always runs on the client, it can be used to present
 652    // the results graphically or save the results to file.
 653
 654    Info("Terminate"," terminate");
 655    AliAnalysisTaskSE::Terminate();
 656
 657    fOutput = dynamic_cast<TList*>(GetOutputData(1));
 658    if (!fOutput) {
 659       printf("ERROR: fOutput not available\n");
 660       return;
 661    }
 662
 663    return;
 664 }
 665
 666 //_______________________________________________________________________________
 667
 668 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
 669 {
 670    //
 671    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 672    //
 673
 674    Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 675
 676    // compute the Delta mass for the D*
 677    if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
 678
 679
 680    fMinMass = mass - range;
 681    fMaxMass = mass + range;
 682
 683    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 684    if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
 685
 686    return;
 687 }
 688
 689 //_______________________________________________________________________________
 690
 691 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
 692 {
 693    //
 694    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 695    //
 696
 697    if (uplimit>lowlimit) {
 698       fMinMass = lowlimit;
 699       fMaxMass = uplimit;
 700    } else {
 701       printf("Error! Lower limit larger than upper limit!\n");
 702       fMinMass = uplimit - uplimit*0.2;
 703       fMaxMass = uplimit;
 704    }
 705
 706    return;
 707 }
 708
 709 //_______________________________________________________________________________
 710
 711 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
 712 {
 713    //
 714    // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
 715    //
 716
 717    if (nptbins>30) {
 718       AliInfo("Maximum number of bins allowed is 30!");
 719       return kFALSE;
 720    }
 721
 722    if (!width) return kFALSE;
 723    for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
 724
 725    return kTRUE;
 726 }
 727
 728 //_______________________________________________________________________________
 729
 730 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
 731 {
 732    //
 733    // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
 734    //
 735
 736    // Statistics
 737    TH1I* hstat = new TH1I("hstat","Statistics",6,-0.5,5.5);
 738    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
 739    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
 740    if(fUseMCInfo){
 741       if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N D");
 742       else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
 743    } else hstat->GetXaxis()->SetBinLabel(3, "N Cand");
 744    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
 745    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
 746    if(fUseMCInfo) {
 747       hstat->GetXaxis()->SetBinLabel(6, "N Background");
 748    }
 749    hstat->SetNdivisions(1);
 750    fOutput->Add(hstat);
 751
 752    TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 753    fOutput->Add(hnCandEv);
 754
 755    // Invariant mass related histograms
 756    const Int_t nbinsmass = 200;
 757    const Int_t ptbinsD=100;
 758    Float_t ptmin=0.,ptmax=50.;
 759    TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
 760    hInvMass->SetStats(kTRUE);
 761    hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
 762    hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 763    fOutput->Add(hInvMass);
 764    if ((fCandidateType==kDstartoKpipi) || fUseMCInfo){
 765       TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 766       fOutput->Add(hnSBCandEv);
 767    }
 768    if (fCandidateType==kDstartoKpipi) {
 769
 770       TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
 771       hPtPion->SetStats(kTRUE);
 772       hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 773       hPtPion->GetYaxis()->SetTitle("entries");
 774       fOutput->Add(hPtPion);
 775    }
 776
 777    fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
 778    fOutput->Add(fhImpPar);
 779
 780    const Int_t nbinsalpha=200;
 781    Float_t minalpha=-TMath::Pi(), maxalpha=TMath::Pi();
 782    const Int_t nbinsdeltaR= 200;
 783    Float_t mindeltaR = 0., maxdeltaR = 10.;
 784    if(fUseMCInfo){
 785       if (fCandidateType==kDstartoKpipi){
 786          TH2F* halphaDDS  =new TH2F("halphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 787          TH2F* halphaDpisS=new TH2F("halphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 788          TH2F* halphaDpiS =new TH2F("halphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 789          TH2F* halphaDKS  =new TH2F("halphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 790
 791          TH2F* halphaDDB  =new TH2F("halphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 792          TH2F* halphaDpisB=new TH2F("halphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 793          TH2F* halphaDpiB =new TH2F("halphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 794          TH2F* halphaDKB  =new TH2F("halphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 795
 796          TH2F* hdeltaRDDS  =new TH2F("hdeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 797          TH2F* hdeltaRDpisS=new TH2F("hdeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 798          TH2F* hdeltaRDpiS =new TH2F("hdeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 799          TH2F* hdeltaRDKS  =new TH2F("hdeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 800
 801          TH2F* hdeltaRDDB  =new TH2F("hdeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 802          TH2F* hdeltaRDpisB=new TH2F("hdeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 803          TH2F* hdeltaRDpiB =new TH2F("hdeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 804          TH2F* hdeltaRDKB  =new TH2F("hdeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 805
 806          fOutput->Add(halphaDDS);
 807          fOutput->Add(halphaDpisS);
 808          fOutput->Add(halphaDpiS);
 809          fOutput->Add(halphaDKS);
 810          fOutput->Add(halphaDDB);
 811          fOutput->Add(halphaDpisB);
 812          fOutput->Add(halphaDpiB);
 813          fOutput->Add(halphaDKB);
 814
 815          fOutput->Add(hdeltaRDDS);
 816          fOutput->Add(hdeltaRDpisS);
 817          fOutput->Add(hdeltaRDpiS);
 818          fOutput->Add(hdeltaRDKS);
 819          fOutput->Add(hdeltaRDDB);
 820          fOutput->Add(hdeltaRDpisB);
 821          fOutput->Add(hdeltaRDpiB);
 822          fOutput->Add(hdeltaRDKB);
 823       }
 824
 825       if (fCandidateType==kD0toKpi){
 826
 827          TH2F* halphaDpiS=new TH2F("halphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 828          TH2F* halphaDKS =new TH2F("halphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 829          TH2F* halphaDpiR=new TH2F("halphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 830          TH2F* halphaDKR =new TH2F("halphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 831
 832          TH2F* halphaDpiB=new TH2F("halphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 833          TH2F* halphaDKB =new TH2F("halphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
 834
 835
 836          TH2F* hdeltaRDpiS=new TH2F("hdeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 837          TH2F* hdeltaRDKS =new TH2F("hdeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 838          TH2F* hdeltaRDpiR=new TH2F("hdeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 839          TH2F* hdeltaRDKR =new TH2F("hdeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 840
 841          TH2F* hdeltaRDpiB=new TH2F("hdeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 842          TH2F* hdeltaRDKB =new TH2F("hdeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
 843
 844          fOutput->Add(halphaDpiS);
 845          fOutput->Add(halphaDKS);
 846          fOutput->Add(halphaDpiR);
 847          fOutput->Add(halphaDKR);
 848          fOutput->Add(halphaDpiB);
 849          fOutput->Add(halphaDKB);
 850
 851          fOutput->Add(hdeltaRDpiS);
 852          fOutput->Add(hdeltaRDKS);
 853          fOutput->Add(hdeltaRDpiR);
 854          fOutput->Add(hdeltaRDKR);
 855          fOutput->Add(hdeltaRDpiB);
 856          fOutput->Add(hdeltaRDKB);
 857
 858       }
 859       fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
 860       fOutput->Add(fhImpParB);
 861
 862       fhInvMassS = new TH1F("hInvMassS", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
 863       fhInvMassS->SetStats(kTRUE);
 864       fhInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
 865       fhInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 866       fOutput->Add(fhInvMassS);
 867
 868       fhInvMassB = new TH1F("hInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
 869       fhInvMassB->SetStats(kTRUE);
 870       fhInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
 871       fhInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 872       fOutput->Add(fhInvMassB);
 873
 874    }
 875    return kTRUE;
 876 }
 877
 878 //_______________________________________________________________________________
 879
 880 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
 881    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
 882
 883    if(!p1 || !p2) return -1;
 884    Double_t phi1=p1->Phi(),eta1=p1->Eta();
 885    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
 886
 887    Double_t dPhi=phi1-phi2;
 888    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
 889    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
 890
 891    Double_t dEta=eta1-eta2;
 892    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
 893    return deltaR;
 894
 895 }