[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)

{
   //
   // 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)
{
   //
   // 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
}

//_______________________________________________________________________________

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; }
   
}

//_______________________________________________________________________________

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
   
   fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
   fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
   
   if (fCandidateType==kDstartoKpipi){
      fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!
      fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
   }
   
   PostData(1, fOutput);
   return;
}

//_______________________________________________________________________________

void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
   //
   // user exec
   //
   
   // add cadidate branch
   fCandidateArray->Delete();
   if (!(InputEvent()->FindListObject(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fCandidateArray);
   if (fCandidateType==kDstartoKpipi){
      fSideBandArray->Delete();
      if (!(InputEvent()->FindListObject(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fSideBandArray);
   }
   //Printf("Arr names %s, %s",fCandidateArray->GetName(),fSideBandArray->GetName());
   // 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;
   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;
   
   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
      charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
      if (!charmCand) continue;
      
      
      if (fUseMCInfo) { // Look in MC, try to simulate the z
      	 if (fCandidateType==kDstartoKpipi) {
      	    AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
      	    mcLabel = temp->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;
      	 }
      	 AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
      	 //if data and Dstar from D0 side band
      	 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*/){	
      	    
      	    new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);
      	    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) {
      	    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){
      	 new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
      	 iSBCand++;
      }
      
      
      Double_t masses[2];
      if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
      	 
      	 //softpion from D* decay
      	 AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
      	 AliAODTrack *track2 = (AliAODTrack*)temp->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 ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {	
      	    masses[0] = temp->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 = temp->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);
   }
   
   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.),
	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	}