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