[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx

/**************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
//
// Task for selecting D mesons to be used as an input for D-Jet correlations
//
//-----------------------------------------------------------------------
// Authors:
// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
// A.Grelli (Utrecht University) a.grelli@uu.nl
// Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
//-----------------------------------------------------------------------

#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TVector3.h>
#include "TROOT.h"
#include <TH3F.h>

#include "AliRDHFCutsDStartoKpipi.h"
#include "AliRDHFCutsD0toKpi.h"
#include "AliAODMCHeader.h"
#include "AliAODHandler.h"
#include "AliAnalysisManager.h"
#include "AliLog.h"
#include "AliAODVertex.h"
#include "AliAODRecoDecay.h"
#include "AliAODRecoCascadeHF.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliESDtrack.h"
#include "AliAODMCParticle.h"
#include "AliAnalysisTaskSEDmesonsFilterCJ.h"

ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)

//_______________________________________________________________________________

AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
AliAnalysisTaskSE(),
fUseMCInfo(kFALSE),
fUseReco(kTRUE),
fCandidateType(0),
fCandidateName(""),
fPDGmother(0),
fNProngs(0),
fBranchName(""),
fOutput(0),
fCuts(0),
fMinMass(0.),
fMaxMass(0.),
fCandidateArray(0),
fSideBandArray(0)

{
   //
   // Default constructor
   //
   
   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
}

//_______________________________________________________________________________

AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
AliAnalysisTaskSE(name),
fUseMCInfo(kFALSE),
fUseReco(kTRUE),
fCandidateType(candtype),
fCandidateName(""),
fPDGmother(0),
fNProngs(0),
fBranchName(""),
fOutput(0),
fCuts(cuts),
fMinMass(0.),
fMaxMass(0.),
fCandidateArray(0),
fSideBandArray(0)
{
   //
   // Constructor. Initialization of Inputs and Outputs
   //
   
   Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
   
   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
   
   const Int_t nptbins = fCuts->GetNPtBins();
   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 };
   
   switch (fCandidateType) {
   case 0 :
      fCandidateName = "D0";
      fPDGmother = 421;
      fNProngs = 2;
      fPDGdaughters[0] = 211;  // pi 
      fPDGdaughters[1] = 321;  // K
      fPDGdaughters[2] = 0;    // empty
      fPDGdaughters[3] = 0;    // empty
      fBranchName = "D0toKpi";
      break;
   case 1 :
      fCandidateName = "Dstar";
      fPDGmother = 413;
      fNProngs = 3;
      fPDGdaughters[1] = 211; // pi soft
      fPDGdaughters[0] = 421; // D0
      fPDGdaughters[2] = 211; // pi fromD0
      fPDGdaughters[3] = 321; // K from D0
      fBranchName = "Dstar";
      
      if (nptbins<=13) {
      	 for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
      } else {
      	 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
      }
      break;
   default :
      printf("%d not accepted!!\n",fCandidateType);
      break;
   }
   
   if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
   if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
   
   DefineOutput(1, TList::Class());       // histos
   DefineOutput(2, AliRDHFCuts::Class()); // my cuts
   DefineOutput(3, TClonesArray::Class()); //array of candidates
   DefineOutput(4, TClonesArray::Class()); //array of SB candidates
}

//_______________________________________________________________________________

AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
{
   //
   // destructor
   //
   
   Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  
   
   if (fOutput) { delete fOutput; fOutput = 0; }
   if (fCuts)   { delete fCuts;   fCuts   = 0; }
   if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
   delete fSideBandArray;
   
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::Init()
{
   //
   // Initialization
   //
   
   if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
   
   switch (fCandidateType) {
   case 0: 
      {
   	 AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
   	 copyfCutsDzero->SetName("AnalysisCutsDzero");
   	 PostData(2, copyfCutsDzero);  // Post the data
      } break;
   case 1: 
      {
      	 AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
      	 copyfCutsDstar->SetName("AnalysisCutsDStar");
      	 PostData(2, copyfCutsDstar); // Post the cuts
      } break;
   default: return;
   }
   
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
{ 
   //
   // output 
   //
   
   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
   
   fOutput = new TList(); fOutput->SetOwner();
   DefineHistoForAnalysis(); // define histograms
   
   if (fCandidateType==kD0toKpi){
      fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
      fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0); 
   }
   
   if (fCandidateType==kDstartoKpipi) {
      fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
      fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); 
   }
   
   fCandidateArray->SetOwner();
   fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
   
   //this is used for the DStar side bands and MC!
   fSideBandArray->SetOwner();
   fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
  
   PostData(1, fOutput);
   PostData(3, fCandidateArray);
   PostData(4, fSideBandArray);
 
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
   //
   // user exec
   //
   // Load the event
   AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
   
   TClonesArray *arrayDStartoD0pi = 0;
   if (!aodEvent && AODEvent() && IsStandardAOD()) {
      
      // In case there is an AOD handler writing a standard AOD, use the AOD 
      // event in memory rather than the input (ESD) event.    
      aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
      
      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
      // have to taken from the AOD event hold by the AliAODExtension
      AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
      if(aodHandler->GetExtensions()) {
      	 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
      	 AliAODEvent *aodFromExt = ext->GetAOD();
      	 arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
      }
   } else {
      arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
   }
   
   if (!arrayDStartoD0pi) {
      AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
      return;
   } else {
      AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
   }
   
   TClonesArray* mcArray = 0x0;
   if (fUseMCInfo) {
      mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
      if (!mcArray) {
      	 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
      	 return;
      }
   }
   
   //Histograms
   TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
   TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
   TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
   TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
   
   TH1F* hPtPion=0x0;
   if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
   hstat->Fill(0);
   
   // fix for temporary bug in ESDfilter 
   // the AODs with null vertex pointer didn't pass the PhysSel
   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
   
   //Event selection
   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
   //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
   if (!iseventselected) return;
   hstat->Fill(1);
   
   const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
   if(fUseReco) hstat->Fill(2, nD);
   
   //D* and D0 prongs needed to MatchToMC method
   Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
   Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
   
   //D0 from D0 bar
   Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K 
   Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi 
   
   Int_t iCand =0;
   Int_t iSBCand=0;
   Int_t isSelected = 0;
   AliAODRecoDecayHF *charmCand = 0;
   AliAODRecoCascadeHF *dstar = 0;
   AliAODMCParticle *charmPart = 0;
   Bool_t isMCBkg=kFALSE;
   
   Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
   
   Int_t mcLabel = -9999;
   Int_t pdgMeson = 413;
   if (fCandidateType==kD0toKpi) pdgMeson = 421;
   
   //clear the TClonesArray from the previous event
   fCandidateArray->Clear();
   fSideBandArray->Clear();
   
   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
      charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
      if (!charmCand) continue;
      
      if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
      
      if (fUseMCInfo) { // Look in MC, try to simulate the z
      	 if (fCandidateType==kDstartoKpipi) {
      	    mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
      	 }
      	 
      	 if (fCandidateType==kD0toKpi) 
      	    mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
      	 
      	 if (mcLabel<=0) isMCBkg=kTRUE;
      	 else hstat->Fill(2);
      	 if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
      }
      
      Double_t ptD = charmCand->Pt();
      
      // region of interest + cuts
      if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;    
      
      if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
      	 //select by track cuts the side band candidates (don't want mass cut)
      	 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent); 
      	 if (!isSelected) continue;
      	 //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
      	 Int_t bin = fCuts->PtBin(ptD);
      	 if(bin<0 || bin>=fCuts->GetNPtBins()) {
      	    AliError(Form("Pt %.3f out of bounds",ptD));
      	    continue;
      	 }
      	 //if data and Dstar from D0 side band
      	 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*/){	
      	    
      	    new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
      	    iSBCand++;
      	 }
      }
      //candidate selected by cuts and PID
      isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
      if (!isSelected) continue;
      
      //for data and MC signal fill fCandidateArray
      if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
      	 // for data or MC with the requirement fUseReco fill with candidates
      	 if(fUseReco) {
      	    if (fCandidateType==kDstartoKpipi){
      	       new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
      	       AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
      	    } else{
      	       new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
      	       //Printf("Filling reco");
      	    }      	    
      	    
      	    hstat->Fill(3);
      	 }
      	 // for MC with requirement particle level fill with AliAODMCParticle
      	 else if (fUseMCInfo) {
      	    new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
      	    //Printf("Filling gen");
      	    hstat->Fill(3);
      	 }
      	 
      	 iCand++;	
      }
      //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
      else if(fUseReco){
      	 if (fCandidateType==kDstartoKpipi){
      	    new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
      	 }
      	 if (fCandidateType==kD0toKpi){
      	    new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
      	 }
      	 iSBCand++;
      }
      
      
      Double_t masses[2];
      if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
      	 
      	 //softpion from D* decay
      	 AliAODTrack *track2 = (AliAODTrack*)dstar->GetBachelor();  
      	 
      	 // select D* in the D0 window.
      	 // In the cut object window is loose to allow for side bands
      	 
      	 
      	 // retrieve the corresponding pt bin for the candidate
      	 // and set the expected D0 width (x3)
      	 //    static const Int_t n = fCuts->GetNPtBins();
      	 Int_t bin = fCuts->PtBin(ptD);
      	 if(bin<0 || bin>=fCuts->GetNPtBins()) {
      	    AliError(Form("Pt %.3f out of bounds",ptD));
      	    continue;
      	 }
      	 
      	 AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
      	 //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
      	 if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {	
      	    masses[0] = dstar->DeltaInvMass(); //D*
      	    masses[1] = 0.; //dummy for D*
      	    
      	    //D*  delta mass
      	    hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
      	    
      	    // D* pt and soft pion pt for good candidates  	      	
      	    Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
      	    Double_t invmassDelta = dstar->DeltaInvMass();
      	    if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
      	 }
      } //Dstar specific
      
      if (fCandidateType==kD0toKpi) { //D0->Kpi
      	 
      	 //needed quantities
      	 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
      	 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
      	 hstat->Fill(3);
      	 
      	 // mass vs pt
      	 if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
      	 if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
      } //D0 specific
      
      charmCand = 0;
      isMCBkg=kFALSE;
   } // end of D cand loop
   
   hnCandEv->Fill(fCandidateArray->GetEntriesFast());
   if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
      Int_t nsbcand=fSideBandArray->GetEntriesFast();
      hstat->Fill(4,nsbcand);
      hnSBCandEv->Fill(nsbcand);
   }
   //Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
   PostData(1, fOutput);
   PostData(3, fCandidateArray);
   PostData(4, fSideBandArray);
  
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
   
   Info("Terminate"," terminate");
   AliAnalysisTaskSE::Terminate();
   
   fOutput = dynamic_cast<TList*>(GetOutputData(1));
   if (!fOutput) {
      printf("ERROR: fOutput not available\n");
      return;
   }
   
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
{
   //
   // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
   //
   
   Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
   
   // compute the Delta mass for the D*
   if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
   
   
   fMinMass = mass - range;
   fMaxMass = mass + range;
   
   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
   if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
   
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
{
   //
   // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
   //
   
   if (uplimit>lowlimit) {
      fMinMass = lowlimit;
      fMaxMass = uplimit;
   } else {
      printf("Error! Lower limit larger than upper limit!\n");
      fMinMass = uplimit - uplimit*0.2;
      fMaxMass = uplimit;
   }
   
   return;
}

//_______________________________________________________________________________

Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
{
   //
   // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
   //
   
   if (nptbins>30) {
      AliInfo("Maximum number of bins allowed is 30!");
      return kFALSE;
   }
   
   if (!width) return kFALSE;
   for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
   
   return kTRUE;
}

//_______________________________________________________________________________

Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
{
   //
   // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
   //
   
   // Statistics 
   TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
   hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
   hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
   if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
   else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
   hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
   if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");  
   hstat->SetNdivisions(1);
   fOutput->Add(hstat);
   
   TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
   fOutput->Add(hnCandEv);
   
   // Invariant mass related histograms
   const Int_t nbinsmass = 200;
   TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
   hInvMass->SetStats(kTRUE);
   hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
   hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
   fOutput->Add(hInvMass);
   
   if (fCandidateType==kDstartoKpipi) {
      TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
      fOutput->Add(hnSBCandEv);
      
      TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
      hPtPion->SetStats(kTRUE);
      hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
      hPtPion->GetYaxis()->SetTitle("entries");
      fOutput->Add(hPtPion);
   }
   
   return kTRUE; 
}
Commit	Line	Data
c683985a	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.),
bbb94467	63	fCandidateArray(0),
bbb94467	64	fSideBandArray(0)
c683985a	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.),
bbb94467	90	fCandidateArray(0),
bbb94467	91	fSideBandArray(0)
c683985a	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
bbb94467	142	DefineOutput(3, TClonesArray::Class()); //array of candidates
bbb94467	143	DefineOutput(4, TClonesArray::Class()); //array of SB candidates
c683985a	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; }
bbb94467	159	delete fSideBandArray;
c683985a	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) {
bbb94467	174	case 0:
bbb94467	175	{
c683985a	176	AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi((static_cast<AliRDHFCutsD0toKpi>(fCuts)));
	177	copyfCutsDzero->SetName("AnalysisCutsDzero");
	178	PostData(2, copyfCutsDzero); // Post the data
bbb94467	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;
c683985a	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
bbb94467	205	if (fCandidateType==kD0toKpi){
	206	fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
	207	fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0);
	208	}
c683985a	209
bbb94467	210	if (fCandidateType==kDstartoKpipi) {
	211	fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
	212	fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0);
c683985a	213	}
c683985a	214
bbb94467	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
c683985a	222	PostData(1, fOutput);
bbb94467	223	PostData(3, fCandidateArray);
	224	PostData(4, fSideBandArray);
	225
c683985a	226	return;
	227	}
	228
	229	//_______________________________________________________________________________
	230
	231	void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
	232	//
	233	// user exec
	234	//
c683985a	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;
bbb94467	308	AliAODRecoCascadeHF *dstar = 0;
c683985a	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
bbb94467	318	//clear the TClonesArray from the previous event
	319	fCandidateArray->Clear();
	320	fSideBandArray->Clear();
	321
c683985a	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
bbb94467	326	if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
c683985a	327
	328	if (fUseMCInfo) { // Look in MC, try to simulate the z
	329	if (fCandidateType==kDstartoKpipi) {
bbb94467	330	mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
c683985a	331	}
	332
	333	if (fCandidateType==kD0toKpi)
	334	mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
	335
	336	if (mcLabel<=0) isMCBkg=kTRUE;
	337	else hstat->Fill(2);
	338	if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
	339	}
	340
	341	Double_t ptD = charmCand->Pt();
	342
	343	// region of interest + cuts
	344	if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
	345
	346	if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
	347	//select by track cuts the side band candidates (don't want mass cut)
	348	isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
	349	if (!isSelected) continue;
	350	//add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
	351	Int_t bin = fCuts->PtBin(ptD);
	352	if(bin<0 \|\| bin>=fCuts->GetNPtBins()) {
	353	AliError(Form("Pt %.3f out of bounds",ptD));
	354	continue;
	355	}
c683985a	356	//if data and Dstar from D0 side band
bbb94467	357	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/){
c683985a	358
bbb94467	359	new ((fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(dstar);
c683985a	360	iSBCand++;
	361	}
	362	}
	363	//candidate selected by cuts and PID
	364	isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
	365	if (!isSelected) continue;
	366
	367	//for data and MC signal fill fCandidateArray
	368	if(!fUseMCInfo \|\| (fUseMCInfo && !isMCBkg)){
	369	// for data or MC with the requirement fUseReco fill with candidates
	370	if(fUseReco) {
bbb94467	371	if (fCandidateType==kDstartoKpipi){
	372	new ((fCandidateArray)[iCand]) AliAODRecoCascadeHF(dstar);
	373	AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
	374	} else{
	375	new ((fCandidateArray)[iCand]) AliAODRecoDecayHF(charmCand);
	376	//Printf("Filling reco");
	377	}
	378
c683985a	379	hstat->Fill(3);
	380	}
	381	// for MC with requirement particle level fill with AliAODMCParticle
	382	else if (fUseMCInfo) {
	383	new ((fCandidateArray)[iCand]) AliAODMCParticle(charmPart);
	384	//Printf("Filling gen");
	385	hstat->Fill(3);
	386	}
	387
	388	iCand++;
	389	}
	390	//for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
	391	else if(fUseReco){
bbb94467	392	if (fCandidateType==kDstartoKpipi){
	393	new ((fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(dstar);
	394	}
	395	if (fCandidateType==kD0toKpi){
	396	new ((fSideBandArray)[iSBCand]) AliAODRecoDecayHF(charmCand);
	397	}
c683985a	398	iSBCand++;
	399	}
	400
	401
	402	Double_t masses[2];
	403	if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
	404
	405	//softpion from D* decay
bbb94467	406	AliAODTrack track2 = (AliAODTrack)dstar->GetBachelor();
c683985a	407
	408	// select D* in the D0 window.
	409	// In the cut object window is loose to allow for side bands
	410
	411
	412	// retrieve the corresponding pt bin for the candidate
	413	// and set the expected D0 width (x3)
	414	// static const Int_t n = fCuts->GetNPtBins();
	415	Int_t bin = fCuts->PtBin(ptD);
	416	if(bin<0 \|\| bin>=fCuts->GetNPtBins()) {
	417	AliError(Form("Pt %.3f out of bounds",ptD));
	418	continue;
	419	}
	420
	421	AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
	422	//consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
bbb94467	423	if ((dstar->InvMassD0()>=(mPDGD0-3.fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.fSigmaD0[bin]))) {
bbb94467	424	masses[0] = dstar->DeltaInvMass(); //D*
c683985a	425	masses[1] = 0.; //dummy for D*
	426
	427	//D* delta mass
	428	hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
	429
	430	// D* pt and soft pion pt for good candidates
	431	Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
bbb94467	432	Double_t invmassDelta = dstar->DeltaInvMass();
c683985a	433	if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
	434	}
	435	} //Dstar specific
	436
	437	if (fCandidateType==kD0toKpi) { //D0->Kpi
	438
	439	//needed quantities
	440	masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
	441	masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
	442	hstat->Fill(3);
	443
	444	// mass vs pt
	445	if (isSelected==1 \|\| isSelected==3) hInvMassptD->Fill(masses[0],ptD);
	446	if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
	447	} //D0 specific
	448
	449	charmCand = 0;
	450	isMCBkg=kFALSE;
	451	} // end of D cand loop
	452
	453	hnCandEv->Fill(fCandidateArray->GetEntriesFast());
	454	if (fCandidateType==kDstartoKpipi \|\| fUseMCInfo) {
	455	Int_t nsbcand=fSideBandArray->GetEntriesFast();
	456	hstat->Fill(4,nsbcand);
	457	hnSBCandEv->Fill(nsbcand);
	458	}
c3c1ab31	459	//Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
bbb94467	460	PostData(1, fOutput);
	461	PostData(3, fCandidateArray);
	462	PostData(4, fSideBandArray);
	463
c683985a	464	return;
	465	}
	466
	467	//_______________________________________________________________________________
	468
	469	void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
	470	{
	471	// The Terminate() function is the last function to be called during
	472	// a query. It always runs on the client, it can be used to present
	473	// the results graphically or save the results to file.
	474
	475	Info("Terminate"," terminate");
	476	AliAnalysisTaskSE::Terminate();
	477
	478	fOutput = dynamic_cast<TList*>(GetOutputData(1));
	479	if (!fOutput) {
	480	printf("ERROR: fOutput not available\n");
	481	return;
	482	}
	483
	484	return;
	485	}
	486
	487	//_______________________________________________________________________________
	488
	489	void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
	490	{
	491	//
	492	// AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
	493	//
	494
	495	Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
	496
	497	// compute the Delta mass for the D*
	498	if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
	499
	500
	501	fMinMass = mass - range;
	502	fMaxMass = mass + range;
	503
	504	AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
	505	if ((fMinMass<0.) \|\| (fMaxMass<=0.) \|\| (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
	506
	507	return;
	508	}
	509
	510	//_______________________________________________________________________________
	511
	512	void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
	513	{
	514	//
	515	// AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
	516	//
	517
	518	if (uplimit>lowlimit) {
	519	fMinMass = lowlimit;
	520	fMaxMass = uplimit;
	521	} else {
	522	printf("Error! Lower limit larger than upper limit!\n");
	523	fMinMass = uplimit - uplimit*0.2;
	524	fMaxMass = uplimit;
	525	}
	526
	527	return;
528	}
529
530	//_______________________________________________________________________________
531
532	Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
533	{
534	//
535	// AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
536	//
537
538	if (nptbins>30) {
539	AliInfo("Maximum number of bins allowed is 30!");
540	return kFALSE;
541	}
542
543	if (!width) return kFALSE;
544	for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
545
546	return kTRUE;
547	}
548
549	//_______________________________________________________________________________
550
551	Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
552	{
553	//
554	// AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
555	//
556
557	// Statistics
558	TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
559	hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
560	hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
561	if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
562	else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
563	hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
564	if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
565	hstat->SetNdivisions(1);
566	fOutput->Add(hstat);
567
568	TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
569	fOutput->Add(hnCandEv);
570
571	// Invariant mass related histograms
572	const Int_t nbinsmass = 200;
573	TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
574	hInvMass->SetStats(kTRUE);
575	hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
576	hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
577	fOutput->Add(hInvMass);
578
579	if (fCandidateType==kDstartoKpipi) {
580	TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
581	fOutput->Add(hnSBCandEv);
582
583	TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
584	hPtPion->SetStats(kTRUE);
585	hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
586	hPtPion->GetYaxis()->SetTitle("entries");
587	fOutput->Add(hPtPion);
588	}
589
590	return kTRUE;
591	}