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
  65 {
  66    //
  67    // Default constructor
  68    //
  69
  70    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
  71    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
  72 }
  73
  74 //_______________________________________________________________________________
  75
  76 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
  77 AliAnalysisTaskSE(name),
  78 fUseMCInfo(kFALSE),
  79 fUseReco(kTRUE),
  80 fCandidateType(candtype),
  81 fCandidateName(""),
  82 fPDGmother(0),
  83 fNProngs(0),
  84 fBranchName(""),
  85 fOutput(0),
  86 fCuts(cuts),
  87 fMinMass(0.),
  88 fMaxMass(0.),
  89 fCandidateArray(0)
  90 {
  91    //
  92    // Constructor. Initialization of Inputs and Outputs
  93    //
  94
  95    Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
  96
  97    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
  98    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
  99
 100    const Int_t nptbins = fCuts->GetNPtBins();
 101    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 };
 102
 103    switch (fCandidateType) {
 104    case 0 :
 105       fCandidateName = "D0";
 106       fPDGmother = 421;
 107       fNProngs = 2;
 108       fPDGdaughters[0] = 211;  // pi
 109       fPDGdaughters[1] = 321;  // K
 110       fPDGdaughters[2] = 0;    // empty
 111       fPDGdaughters[3] = 0;    // empty
 112       fBranchName = "D0toKpi";
 113       break;
 114    case 1 :
 115       fCandidateName = "Dstar";
 116       fPDGmother = 413;
 117       fNProngs = 3;
 118       fPDGdaughters[1] = 211; // pi soft
 119       fPDGdaughters[0] = 421; // D0
 120       fPDGdaughters[2] = 211; // pi fromD0
 121       fPDGdaughters[3] = 321; // K from D0
 122       fBranchName = "Dstar";
 123
 124       if (nptbins<=13) {
 125          for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
 126       } else {
 127          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 128       }
 129       break;
 130    default :
 131       printf("%d not accepted!!\n",fCandidateType);
 132       break;
 133    }
 134
 135    if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
 136    if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 137
 138    DefineOutput(1, TList::Class());       // histos
 139    DefineOutput(2, AliRDHFCuts::Class()); // my cuts
 140 }
 141
 142 //_______________________________________________________________________________
 143
 144 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
 145 {
 146    //
 147    // destructor
 148    //
 149
 150    Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
 151
 152    if (fOutput) { delete fOutput; fOutput = 0; }
 153    if (fCuts)   { delete fCuts;   fCuts   = 0; }
 154    if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
 155
 156 }
 157
 158 //_______________________________________________________________________________
 159
 160 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
 161 {
 162    //
 163    // Initialization
 164    //
 165
 166    if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
 167
 168    switch (fCandidateType) {
 169    case 0: {
 170          AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 171          copyfCutsDzero->SetName("AnalysisCutsDzero");
 172          PostData(2, copyfCutsDzero);  // Post the data
 173    } break;
 174 case 1: {
 175       AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 176       copyfCutsDstar->SetName("AnalysisCutsDStar");
 177       PostData(2, copyfCutsDstar); // Post the cuts
 178 } break;
 179 default: return;
 180    }
 181
 182    return;
 183 }
 184
 185 //_______________________________________________________________________________
 186
 187 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
 188 {
 189    //
 190    // output
 191    //
 192
 193    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 194
 195    fOutput = new TList(); fOutput->SetOwner();
 196    DefineHistoForAnalysis(); // define histograms
 197
 198    fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
 199    fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
 200
 201    if (fCandidateType==kDstartoKpipi){
 202       fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!
 203       fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
 204    }
 205
 206    PostData(1, fOutput);
 207    return;
 208 }
 209
 210 //_______________________________________________________________________________
 211
 212 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
 213    //
 214    // user exec
 215    //
 216
 217    // add cadidate branch
 218    fCandidateArray->Delete();
 219    if (!(InputEvent()->FindListObject(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fCandidateArray);
 220    if (fCandidateType==kDstartoKpipi){
 221       fSideBandArray->Delete();
 222       if (!(InputEvent()->FindListObject(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fSideBandArray);
 223    }
 224    //Printf("Arr names %s, %s",fCandidateArray->GetName(),fSideBandArray->GetName());
 225    // Load the event
 226    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 227
 228    TClonesArray *arrayDStartoD0pi = 0;
 229    if (!aodEvent && AODEvent() && IsStandardAOD()) {
 230
 231       // In case there is an AOD handler writing a standard AOD, use the AOD
 232       // event in memory rather than the input (ESD) event.
 233       aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
 234
 235       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 236       // have to taken from the AOD event hold by the AliAODExtension
 237       AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 238       if(aodHandler->GetExtensions()) {
 239          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 240          AliAODEvent *aodFromExt = ext->GetAOD();
 241          arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 242       }
 243    } else {
 244       arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 245    }
 246
 247    if (!arrayDStartoD0pi) {
 248       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 249       return;
 250    } else {
 251       AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 252    }
 253
 254    TClonesArray* mcArray = 0x0;
 255    if (fUseMCInfo) {
 256       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 257       if (!mcArray) {
 258          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 259          return;
 260       }
 261    }
 262
 263    //Histograms
 264    TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
 265    TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
 266    TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
 267    TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
 268
 269    TH1F* hPtPion=0x0;
 270    if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
 271    hstat->Fill(0);
 272
 273    // fix for temporary bug in ESDfilter
 274    // the AODs with null vertex pointer didn't pass the PhysSel
 275    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 276
 277    //Event selection
 278    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 279    //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 280    if (!iseventselected) return;
 281    hstat->Fill(1);
 282
 283    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
 284    if(fUseReco) hstat->Fill(2, nD);
 285
 286    //D* and D0 prongs needed to MatchToMC method
 287    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
 288    Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
 289
 290    //D0 from D0 bar
 291    Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K
 292    Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi
 293
 294    Int_t iCand =0;
 295    Int_t iSBCand=0;
 296    Int_t isSelected = 0;
 297    AliAODRecoDecayHF *charmCand = 0;
 298    AliAODMCParticle *charmPart = 0;
 299    Bool_t isMCBkg=kFALSE;
 300
 301    Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 302
 303    Int_t mcLabel = -9999;
 304    Int_t pdgMeson = 413;
 305    if (fCandidateType==kD0toKpi) pdgMeson = 421;
 306
 307    for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
 308       charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
 309       if (!charmCand) continue;
 310
 311
 312       if (fUseMCInfo) { // Look in MC, try to simulate the z
 313          if (fCandidateType==kDstartoKpipi) {
 314             AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 315             mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 316          }
 317
 318          if (fCandidateType==kD0toKpi)
 319             mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
 320
 321          if (mcLabel<=0) isMCBkg=kTRUE;
 322          else hstat->Fill(2);
 323          if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
 324       }
 325
 326       Double_t ptD = charmCand->Pt();
 327
 328       // region of interest + cuts
 329       if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
 330
 331       if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
 332          //select by track cuts the side band candidates (don't want mass cut)
 333          isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
 334          if (!isSelected) continue;
 335          //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
 336          Int_t bin = fCuts->PtBin(ptD);
 337          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 338             AliError(Form("Pt %.3f out of bounds",ptD));
 339             continue;
 340          }
 341          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 342          //if data and Dstar from D0 side band
 343          if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
 344
 345             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);
 346             iSBCand++;
 347          }
 348       }
 349       //candidate selected by cuts and PID
 350       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
 351       if (!isSelected) continue;
 352
 353       //for data and MC signal fill fCandidateArray
 354       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
 355          // for data or MC with the requirement fUseReco fill with candidates
 356          if(fUseReco) {
 357             new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
 358             //Printf("Filling reco");
 359             hstat->Fill(3);
 360          }
 361          // for MC with requirement particle level fill with AliAODMCParticle
 362          else if (fUseMCInfo) {
 363             new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
 364             //Printf("Filling gen");
 365             hstat->Fill(3);
 366          }
 367
 368          iCand++;
 369       }
 370       //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
 371       else if(fUseReco){
 372          new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
 373          iSBCand++;
 374       }
 375
 376
 377       Double_t masses[2];
 378       if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
 379
 380          //softpion from D* decay
 381          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 382          AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();
 383
 384          // select D* in the D0 window.
 385          // In the cut object window is loose to allow for side bands
 386
 387
 388          // retrieve the corresponding pt bin for the candidate
 389          // and set the expected D0 width (x3)
 390          //    static const Int_t n = fCuts->GetNPtBins();
 391          Int_t bin = fCuts->PtBin(ptD);
 392          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 393             AliError(Form("Pt %.3f out of bounds",ptD));
 394             continue;
 395          }
 396
 397          AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
 398          //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
 399          if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
 400             masses[0] = temp->DeltaInvMass(); //D*
 401             masses[1] = 0.; //dummy for D*
 402
 403             //D*  delta mass
 404             hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
 405
 406             // D* pt and soft pion pt for good candidates
 407             Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 408             Double_t invmassDelta = temp->DeltaInvMass();
 409             if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
 410          }
 411       } //Dstar specific
 412
 413       if (fCandidateType==kD0toKpi) { //D0->Kpi
 414
 415          //needed quantities
 416          masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 417          masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 418          hstat->Fill(3);
 419
 420          // mass vs pt
 421          if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
 422          if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
 423       } //D0 specific
 424
 425       charmCand = 0;
 426       isMCBkg=kFALSE;
 427    } // end of D cand loop
 428
 429    hnCandEv->Fill(fCandidateArray->GetEntriesFast());
 430    if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
 431       Int_t nsbcand=fSideBandArray->GetEntriesFast();
 432       hstat->Fill(4,nsbcand);
 433       hnSBCandEv->Fill(nsbcand);
 434    }
 435
 436    return;
 437 }
 438
 439 //_______________________________________________________________________________
 440
 441 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
 442 {
 443    // The Terminate() function is the last function to be called during
 444    // a query. It always runs on the client, it can be used to present
 445    // the results graphically or save the results to file.
 446
 447    Info("Terminate"," terminate");
 448    AliAnalysisTaskSE::Terminate();
 449
 450    fOutput = dynamic_cast<TList*>(GetOutputData(1));
 451    if (!fOutput) {
 452       printf("ERROR: fOutput not available\n");
 453       return;
 454    }
 455
 456    return;
 457 }
 458
 459 //_______________________________________________________________________________
 460
 461 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
 462 {
 463    //
 464    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 465    //
 466
 467    Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 468
 469    // compute the Delta mass for the D*
 470    if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
 471
 472
 473    fMinMass = mass - range;
 474    fMaxMass = mass + range;
 475
 476    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 477    if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
 478
 479    return;
 480 }
 481
 482 //_______________________________________________________________________________
 483
 484 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
 485 {
 486    //
 487    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 488    //
 489
 490    if (uplimit>lowlimit) {
 491       fMinMass = lowlimit;
 492       fMaxMass = uplimit;
 493    } else {
 494       printf("Error! Lower limit larger than upper limit!\n");
 495       fMinMass = uplimit - uplimit*0.2;
 496       fMaxMass = uplimit;
 497    }
 498
 499    return;
 500 }
 501
 502 //_______________________________________________________________________________
 503
 504 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
 505 {
 506    //
 507    // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
 508    //
 509
 510    if (nptbins>30) {
 511       AliInfo("Maximum number of bins allowed is 30!");
 512       return kFALSE;
 513    }
 514
 515    if (!width) return kFALSE;
 516    for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
 517
 518    return kTRUE;
 519 }
 520
 521 //_______________________________________________________________________________
 522
 523 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
 524 {
 525    //
 526    // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
 527    //
 528
 529    // Statistics
 530    TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
 531    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
 532    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
 533    if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
 534    else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
 535    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
 536    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
 537    hstat->SetNdivisions(1);
 538    fOutput->Add(hstat);
 539
 540    TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 541    fOutput->Add(hnCandEv);
 542
 543    // Invariant mass related histograms
 544    const Int_t nbinsmass = 200;
 545    TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
 546    hInvMass->SetStats(kTRUE);
 547    hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
 548    hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 549    fOutput->Add(hInvMass);
 550
 551    if (fCandidateType==kDstartoKpipi) {
 552       TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 553       fOutput->Add(hnSBCandEv);
 554
 555       TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
 556       hPtPion->SetStats(kTRUE);
 557       hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 558       hPtPion->GetYaxis()->SetTitle("entries");
 559       fOutput->Add(hPtPion);
 560    }
 561
 562    return kTRUE;
 563 }