[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDs.cxx

/**************************************************************************
 * Copyright(c) 2007-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.                  *
 **************************************************************************/

/* $Id: AliITSCorrMapSDD.cxx 32906 2009-06-12 16:56:53Z prino $ */

///////////////////////////////////////////////////////////////////
//                                                               //
//  Analysis task to produce Ds candidates mass spectra          //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TNtuple.h>
#include <TList.h>
#include <TString.h>
#include <TH1F.h>
#include <TMath.h>
#include <TDatabasePDG.h>

#include "AliAnalysisManager.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODRecoDecayHF3Prong.h"
#include "AliAnalysisVertexingHF.h"
#include "AliRDHFCutsDstoKKpi.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSEDs.h"

ClassImp(AliAnalysisTaskSEDs)


//________________________________________________________________________
AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
  AliAnalysisTaskSE(),
  fOutput(0), 
  fHistNEvents(0),
  fReadMC(kFALSE),
  fNPtBins(0),
  fListCuts(0),
  fMassRange(0.2),
  fMassBinSize(0.002),
  fProdCuts(0),
  fAnalysisCuts(0)
{
  // Default constructor
}

//________________________________________________________________________
AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts):
  AliAnalysisTaskSE(name),
  fOutput(0),
  fHistNEvents(0),
  fReadMC(kFALSE),
  fNPtBins(0),
  fListCuts(0),
  fMassRange(0.2),
  fMassBinSize(0.002),
  fProdCuts(productioncuts),
  fAnalysisCuts(analysiscuts)
{
  // Default constructor
  // Output slot #1 writes into a TList container
  Int_t nptbins=fAnalysisCuts->GetNPtBins();
  Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
  SetPtBins(nptbins,ptlim);
 
  DefineOutput(1,TList::Class());  //My private output

  DefineOutput(2,TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
  // define pt bins for analysis
  if(n>kMaxPtBins){
    printf("Max. number of Pt bins = %d\n",kMaxPtBins);
    fNPtBins=kMaxPtBins;
    fPtLimits[0]=0.;
    fPtLimits[1]=1.;
    fPtLimits[2]=3.;
    fPtLimits[3]=5.;
    fPtLimits[4]=10.;
    for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
  }else{
    fNPtBins=n;
    for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
    for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
  }
  if(fDebug > 1){
    printf("Number of Pt bins = %d\n",fNPtBins);
    for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);    
  }
}
//________________________________________________________________________
AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
{
  // Destructor
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
  if (fListCuts) {
    delete fListCuts;
    fListCuts = 0;
  }

  if (fProdCuts) {
    delete fProdCuts;
    fProdCuts = 0;
  }
  if (fAnalysisCuts) {
    delete fAnalysisCuts;
    fAnalysisCuts = 0;
  }
}  

//________________________________________________________________________
void AliAnalysisTaskSEDs::Init()
{
  // Initialization

  if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");

  fListCuts=new TList();
  AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi();
  production=fProdCuts;
  AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi();
  analysis=fAnalysisCuts;
  
  fListCuts->Add(production);
  fListCuts->Add(analysis);
  PostData(2,fListCuts);
  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDs::UserCreateOutputObjects()
{
  // Create the output container
  //
  if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList();
  fOutput->SetOwner();
  fOutput->SetName("OutputHistos");

  Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
  Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
  if(nInvMassBins%2==1) nInvMassBins++;
  Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
  Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;

  TString hisname;
  Int_t index;
  for(Int_t i=0;i<fNPtBins;i++){
    index=GetHistoIndex(i);
    hisname.Form("hMassAllPt%d",i);
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHist[index]->Sumw2();
    hisname.Form("hMassAllPt%dCuts",i);
    fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHistCuts[index]->Sumw2();
    hisname.Form("hCosPAllPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenAllPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hDalitzAllPt%d",i);
    fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
    fDalitz[index]->Sumw2();

    index=GetSignalHistoIndex(i);    
    hisname.Form("hMassSigPt%d",i);
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHist[index]->Sumw2();
    hisname.Form("hMassSigPt%dCuts",i);
    fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHistCuts[index]->Sumw2();
    hisname.Form("hCosPSigPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenSigPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hDalitzSigPt%d",i);
    fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
    fDalitz[index]->Sumw2();

    hisname.Form("hMassBkgPt%d",i);
    index=GetBackgroundHistoIndex(i);    
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHist[index]->Sumw2();
    hisname.Form("hMassBkgPt%dCuts",i);
    fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
    fMassHistCuts[index]->Sumw2();
    hisname.Form("hCosPBkgPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenBkgPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hDalitzBkgPt%d",i);
    fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
    fDalitz[index]->Sumw2();
  }

  for(Int_t i=0; i<3*fNPtBins; i++){
    fOutput->Add(fMassHist[i]);
    fOutput->Add(fMassHistCuts[i]);
    fOutput->Add(fCosPHist[i]);
    fOutput->Add(fDLenHist[i]);
    fOutput->Add(fDalitz[i]);
  }

  fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
  fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
  fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
  fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
  fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
  fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
  for(Int_t i=0;i<3;i++){
    fChanHist[i]->Sumw2();
    fChanHist[i]->SetMinimum(0);
    fChanHistCuts[i]->Sumw2();
    fChanHistCuts[i]->SetMinimum(0);
    fOutput->Add(fChanHist[i]);
    fOutput->Add(fChanHistCuts[i]);
  }

  fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
  fHistNEvents->Sumw2();
  fHistNEvents->SetMinimum(0);
  fOutput->Add(fHistNEvents);


  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
{
  // Ds selection for current event, fill mass histos and selecetion variable histo
  // separate signal and backgound if fReadMC is activated

  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
  fHistNEvents->Fill(0); // count event
  // Post the data already here
  PostData(1,fOutput);
  

  TClonesArray *array3Prong = 0;
  if(!aod && 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.    
    aod = 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();
      array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
    }
  } else {
    array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
  }

  if(!array3Prong) {
    printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
    return;
  }

 
  TClonesArray *arrayMC=0;
  AliAODMCHeader *mcHeader=0;

  // AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  //    vtx1->Print();
  
  // load MC particles
  if(fReadMC){
    
    arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
    if(!arrayMC) {
      printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
      return;
    }
    
  // load MC header
    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
    if(!mcHeader) {
      printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
      return;
    }
  }
  
  Int_t n3Prong = array3Prong->GetEntriesFast();
  if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
  
  
  Int_t pdgDstoKKpi[3]={321,321,211};
  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
    
    
    Bool_t unsetvtx=kFALSE;
    if(!d->GetOwnPrimaryVtx()){
      d->SetOwnPrimaryVtx(vtx1);
      unsetvtx=kTRUE;
    }

    Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
    if(retCodeProductionCuts>0){
      Int_t isKKpi=retCodeProductionCuts&1;
      Int_t ispiKK=retCodeProductionCuts&2;
//       Int_t isPhi=retCodeProductionCuts&4;
//       Int_t isK0star=retCodeProductionCuts&8;
      Double_t ptCand = d->Pt();
      Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
      Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
      Int_t isKKpiAC=retCodeAnalysisCuts&1;
      Int_t ispiKKAC=retCodeAnalysisCuts&2;
//       Int_t isPhiAC=retCodeAnalysisCuts&4;
//       Int_t isK0starAC=retCodeAnalysisCuts&8;

      Int_t labDs=-1;
      if(fReadMC){
	labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
      }

      Double_t dlen=d->DecayLength();
      Double_t cosp=d->CosPointingAngle();
      Int_t index=GetHistoIndex(iPtBin);
      Int_t type=0;
      if(isKKpi) type+=1;
      if(ispiKK) type+=2;
      Int_t typeAC=0;
      if(isKKpiAC) typeAC+=1;
      if(ispiKKAC) typeAC+=2;
      fCosPHist[index]->Fill(cosp);
      fDLenHist[index]->Fill(dlen);
      fChanHist[0]->Fill(type);
      if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
      if(fReadMC){
	if(labDs>=0) {	  
	  index=GetSignalHistoIndex(iPtBin);
	  fCosPHist[index]->Fill(cosp);
	  fDLenHist[index]->Fill(dlen);
	  fChanHist[1]->Fill(type);	  
	  if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
	}else{
	  index=GetBackgroundHistoIndex(iPtBin);
	  fCosPHist[index]->Fill(cosp);
	  fDLenHist[index]->Fill(dlen);
	  fChanHist[2]->Fill(type);	  
	  if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
	}
      }
      if(isKKpi){
	index=GetHistoIndex(iPtBin);
	Double_t invMass=d->InvMassDsKKpi();
	fMassHist[index]->Fill(invMass);
	Double_t mass01=d->InvMass2Prongs(0,1,321,321);
	Double_t mass12=d->InvMass2Prongs(1,2,321,211);
	fDalitz[index]->Fill(mass01,mass12);
	if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
	if(fReadMC){
	  if(labDs>=0) {	  
	    index=GetSignalHistoIndex(iPtBin);
	    fMassHist[index]->Fill(invMass);
	    fDalitz[index]->Fill(mass01,mass12);
	    if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
	  }else{
	    index=GetBackgroundHistoIndex(iPtBin);
	    fMassHist[index]->Fill(invMass);
	    fDalitz[index]->Fill(mass01,mass12);
	    if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
	  }
	}	
      }
      if(ispiKK){
	index=GetHistoIndex(iPtBin);
	Double_t invMass=d->InvMassDspiKK();
	fMassHist[index]->Fill(invMass);
	if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
	if(fReadMC){
	  if(labDs>=0) {	  
	    index=GetSignalHistoIndex(iPtBin);
	    fMassHist[index]->Fill(invMass);
	    if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
	  }else{
	    index=GetBackgroundHistoIndex(iPtBin);
	    fMassHist[index]->Fill(invMass);
	    if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
	  }
	}
      }
    }
    if(unsetvtx) d->UnsetOwnPrimaryVtx();
  }
 
   
  PostData(1,fOutput);    
  return;
}

//_________________________________________________________________

void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    printf("ERROR: fOutput not available\n");
    return;
  }
  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
  fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
  fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
  fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
  fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
  fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
  fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));


  TString hisname;
  Int_t index;
  for(Int_t i=0;i<fNPtBins;i++){

    index=GetHistoIndex(i);
    hisname.Form("hMassAllPt%d",i);
    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hMassAllPt%dCuts",i);
    fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hCosPAllPt%d",i);
    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDLenAllPt%d",i);
    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDalitzAllPt%d",i);
    fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));

    index=GetSignalHistoIndex(i);    
    hisname.Form("hMassSigPt%d",i);
    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hMassSigPt%dCuts",i);
    fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hCosPSigPt%d",i);
    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDLenSigPt%d",i);
    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDalitzSigPt%d",i);
    fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));

    index=GetBackgroundHistoIndex(i);    
    hisname.Form("hMassBkgPt%d",i);
    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hMassBkgPt%dCuts",i);
    fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hCosPBkgPt%d",i);
    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDLenBkgPt%d",i);
    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
    hisname.Form("hDalitzBkgPt%d",i);
    fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));

  }
  return;
}
Commit	Line	Data
25c94fc8	1	/**************************************************************************
	2	* Copyright(c) 2007-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	/* $Id: AliITSCorrMapSDD.cxx 32906 2009-06-12 16:56:53Z prino $ */
	17
	18	///////////////////////////////////////////////////////////////////
	19	// //
	20	// Analysis task to produce Ds candidates mass spectra //
	21	// Origin: F.Prino, Torino, prino@to.infn.it //
	22	// //
	23	///////////////////////////////////////////////////////////////////
	24
	25	#include <TClonesArray.h>
	26	#include <TNtuple.h>
	27	#include <TList.h>
	28	#include <TString.h>
	29	#include <TH1F.h>
	30	#include <TMath.h>
	31	#include <TDatabasePDG.h>
	32
	33	#include "AliAnalysisManager.h"
	34	#include "AliAODHandler.h"
	35	#include "AliAODEvent.h"
	36	#include "AliAODVertex.h"
	37	#include "AliAODTrack.h"
	38	#include "AliAODMCHeader.h"
	39	#include "AliAODMCParticle.h"
	40	#include "AliAODRecoDecayHF3Prong.h"
	41	#include "AliAnalysisVertexingHF.h"
4c230f6d	42	#include "AliRDHFCutsDstoKKpi.h"
25c94fc8	43	#include "AliAnalysisTaskSE.h"
	44	#include "AliAnalysisTaskSEDs.h"
	45
	46	ClassImp(AliAnalysisTaskSEDs)
	47
	48
	49	//________________________________________________________________________
	50	AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
	51	AliAnalysisTaskSE(),
	52	fOutput(0),
	53	fHistNEvents(0),
	54	fReadMC(kFALSE),
	55	fNPtBins(0),
4c230f6d	56	fListCuts(0),
25c94fc8	57	fMassRange(0.2),
a0bb2d06	58	fMassBinSize(0.002),
4c230f6d	59	fProdCuts(0),
4c230f6d	60	fAnalysisCuts(0)
25c94fc8	61	{
	62	// Default constructor
	63	}
	64
	65	//________________________________________________________________________
4c230f6d	66	AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char name, AliRDHFCutsDstoKKpi productioncuts, AliRDHFCutsDstoKKpi* analysiscuts):
25c94fc8	67	AliAnalysisTaskSE(name),
	68	fOutput(0),
	69	fHistNEvents(0),
	70	fReadMC(kFALSE),
	71	fNPtBins(0),
4c230f6d	72	fListCuts(0),
25c94fc8	73	fMassRange(0.2),
a0bb2d06	74	fMassBinSize(0.002),
4c230f6d	75	fProdCuts(productioncuts),
4c230f6d	76	fAnalysisCuts(analysiscuts)
25c94fc8	77	{
	78	// Default constructor
	79	// Output slot #1 writes into a TList container
4c230f6d	80	Int_t nptbins=fAnalysisCuts->GetNPtBins();
	81	Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
	82	SetPtBins(nptbins,ptlim);
25c94fc8	83
	84	DefineOutput(1,TList::Class()); //My private output
	85
4c230f6d	86	DefineOutput(2,TList::Class());
25c94fc8	87	}
	88
	89	//________________________________________________________________________
4c230f6d	90	void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
25c94fc8	91	// define pt bins for analysis
	92	if(n>kMaxPtBins){
	93	printf("Max. number of Pt bins = %d\n",kMaxPtBins);
	94	fNPtBins=kMaxPtBins;
	95	fPtLimits[0]=0.;
	96	fPtLimits[1]=1.;
	97	fPtLimits[2]=3.;
	98	fPtLimits[3]=5.;
	99	fPtLimits[4]=10.;
	100	for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
	101	}else{
	102	fNPtBins=n;
	103	for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
	104	for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
	105	}
	106	if(fDebug > 1){
	107	printf("Number of Pt bins = %d\n",fNPtBins);
	108	for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
	109	}
	110	}
	111	//________________________________________________________________________
	112	AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
	113	{
	114	// Destructor
	115	if (fOutput) {
	116	delete fOutput;
	117	fOutput = 0;
	118	}
4c230f6d	119	if (fListCuts) {
	120	delete fListCuts;
	121	fListCuts = 0;
	122	}
	123
	124	if (fProdCuts) {
	125	delete fProdCuts;
	126	fProdCuts = 0;
	127	}
	128	if (fAnalysisCuts) {
	129	delete fAnalysisCuts;
	130	fAnalysisCuts = 0;
25c94fc8	131	}
	132	}
	133
	134	//________________________________________________________________________
	135	void AliAnalysisTaskSEDs::Init()
	136	{
	137	// Initialization
	138
	139	if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
	140
4c230f6d	141	fListCuts=new TList();
	142	AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi();
	143	production=fProdCuts;
	144	AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi();
	145	analysis=fAnalysisCuts;
	146
	147	fListCuts->Add(production);
	148	fListCuts->Add(analysis);
	149	PostData(2,fListCuts);
25c94fc8	150	return;
	151	}
	152
	153	//________________________________________________________________________
	154	void AliAnalysisTaskSEDs::UserCreateOutputObjects()
	155	{
	156	// Create the output container
	157	//
	158	if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
	159
	160	// Several histograms are more conveniently managed in a TList
	161	fOutput = new TList();
	162	fOutput->SetOwner();
	163	fOutput->SetName("OutputHistos");
	164
	165	Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
a0bb2d06	166	Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
	167	if(nInvMassBins%2==1) nInvMassBins++;
	168	Double_t minMass=massDs-0.5nInvMassBinsfMassBinSize;
	169	Double_t maxMass=massDs+0.5nInvMassBinsfMassBinSize;
	170
25c94fc8	171	TString hisname;
	172	Int_t index;
	173	for(Int_t i=0;i<fNPtBins;i++){
	174	index=GetHistoIndex(i);
	175	hisname.Form("hMassAllPt%d",i);
a0bb2d06	176	fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	177	fMassHist[index]->Sumw2();
25c94fc8	178	hisname.Form("hMassAllPt%dCuts",i);
a0bb2d06	179	fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	180	fMassHistCuts[index]->Sumw2();
	181	hisname.Form("hCosPAllPt%d",i);
	182	fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
	183	fCosPHist[index]->Sumw2();
	184	hisname.Form("hDLenAllPt%d",i);
	185	fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
	186	fDLenHist[index]->Sumw2();
	187	hisname.Form("hDalitzAllPt%d",i);
	188	fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
	189	fDalitz[index]->Sumw2();
	190
	191	index=GetSignalHistoIndex(i);
	192	hisname.Form("hMassSigPt%d",i);
a0bb2d06	193	fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	194	fMassHist[index]->Sumw2();
25c94fc8	195	hisname.Form("hMassSigPt%dCuts",i);
a0bb2d06	196	fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	197	fMassHistCuts[index]->Sumw2();
	198	hisname.Form("hCosPSigPt%d",i);
	199	fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
	200	fCosPHist[index]->Sumw2();
	201	hisname.Form("hDLenSigPt%d",i);
	202	fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
	203	fDLenHist[index]->Sumw2();
	204	hisname.Form("hDalitzSigPt%d",i);
	205	fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
	206	fDalitz[index]->Sumw2();
	207
	208	hisname.Form("hMassBkgPt%d",i);
	209	index=GetBackgroundHistoIndex(i);
a0bb2d06	210	fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	211	fMassHist[index]->Sumw2();
25c94fc8	212	hisname.Form("hMassBkgPt%dCuts",i);
a0bb2d06	213	fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
25c94fc8	214	fMassHistCuts[index]->Sumw2();
	215	hisname.Form("hCosPBkgPt%d",i);
	216	fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
	217	fCosPHist[index]->Sumw2();
	218	hisname.Form("hDLenBkgPt%d",i);
	219	fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
	220	fDLenHist[index]->Sumw2();
	221	hisname.Form("hDalitzBkgPt%d",i);
	222	fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
	223	fDalitz[index]->Sumw2();
	224	}
	225
	226	for(Int_t i=0; i<3*fNPtBins; i++){
	227	fOutput->Add(fMassHist[i]);
	228	fOutput->Add(fMassHistCuts[i]);
	229	fOutput->Add(fCosPHist[i]);
	230	fOutput->Add(fDLenHist[i]);
	231	fOutput->Add(fDalitz[i]);
	232	}
	233
	234	fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
	235	fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
	236	fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
	237	fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
	238	fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
	239	fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
	240	for(Int_t i=0;i<3;i++){
	241	fChanHist[i]->Sumw2();
	242	fChanHist[i]->SetMinimum(0);
	243	fChanHistCuts[i]->Sumw2();
	244	fChanHistCuts[i]->SetMinimum(0);
	245	fOutput->Add(fChanHist[i]);
	246	fOutput->Add(fChanHistCuts[i]);
	247	}
	248
	249	fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
	250	fHistNEvents->Sumw2();
	251	fHistNEvents->SetMinimum(0);
	252	fOutput->Add(fHistNEvents);
	253
	254
	255	return;
	256	}
	257
	258	//________________________________________________________________________
	259	void AliAnalysisTaskSEDs::UserExec(Option_t /option*/)
	260	{
	261	// Ds selection for current event, fill mass histos and selecetion variable histo
	262	// separate signal and backgound if fReadMC is activated
	263
	264	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	265	fHistNEvents->Fill(0); // count event
	266	// Post the data already here
	267	PostData(1,fOutput);
	268
	269
	270	TClonesArray *array3Prong = 0;
	271	if(!aod && AODEvent() && IsStandardAOD()) {
	272	// In case there is an AOD handler writing a standard AOD, use the AOD
	273	// event in memory rather than the input (ESD) event.
	274	aod = dynamic_cast<AliAODEvent*> (AODEvent());
	275	// in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
	276	// have to taken from the AOD event hold by the AliAODExtension
	277	AliAODHandler* aodHandler = (AliAODHandler*)
278	((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
279	if(aodHandler->GetExtensions()) {
280	AliAODExtension ext = (AliAODExtension)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
281	AliAODEvent *aodFromExt = ext->GetAOD();
282	array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
283	}
284	} else {
285	array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
286	}
287
288	if(!array3Prong) {
289	printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
290	return;
291	}
292
293
294	TClonesArray *arrayMC=0;
295	AliAODMCHeader *mcHeader=0;
296
297	// AOD primary vertex
298	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
299	// vtx1->Print();
300
301	// load MC particles
302	if(fReadMC){
303
304	arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
305	if(!arrayMC) {
306	printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
307	return;
308	}
309
310	// load MC header
311	mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
312	if(!mcHeader) {
313	printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
314	return;
315	}
316	}
317
318	Int_t n3Prong = array3Prong->GetEntriesFast();
319	if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
320
321
322	Int_t pdgDstoKKpi[3]={321,321,211};
25c94fc8	323	for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
	324	AliAODRecoDecayHF3Prong d = (AliAODRecoDecayHF3Prong)array3Prong->UncheckedAt(i3Prong);
	325
	326
	327	Bool_t unsetvtx=kFALSE;
	328	if(!d->GetOwnPrimaryVtx()){
	329	d->SetOwnPrimaryVtx(vtx1);
	330	unsetvtx=kTRUE;
	331	}
4c230f6d	332
	333	Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
	334	if(retCodeProductionCuts>0){
	335	Int_t isKKpi=retCodeProductionCuts&1;
	336	Int_t ispiKK=retCodeProductionCuts&2;
	337	// Int_t isPhi=retCodeProductionCuts&4;
	338	// Int_t isK0star=retCodeProductionCuts&8;
25c94fc8	339	Double_t ptCand = d->Pt();
4c230f6d	340	Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
	341	Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
	342	Int_t isKKpiAC=retCodeAnalysisCuts&1;
	343	Int_t ispiKKAC=retCodeAnalysisCuts&2;
	344	// Int_t isPhiAC=retCodeAnalysisCuts&4;
	345	// Int_t isK0starAC=retCodeAnalysisCuts&8;
	346
25c94fc8	347	Int_t labDs=-1;
	348	if(fReadMC){
	349	labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
	350	}
	351
	352	Double_t dlen=d->DecayLength();
	353	Double_t cosp=d->CosPointingAngle();
	354	Int_t index=GetHistoIndex(iPtBin);
	355	Int_t type=0;
	356	if(isKKpi) type+=1;
	357	if(ispiKK) type+=2;
4c230f6d	358	Int_t typeAC=0;
	359	if(isKKpiAC) typeAC+=1;
	360	if(ispiKKAC) typeAC+=2;
25c94fc8	361	fCosPHist[index]->Fill(cosp);
	362	fDLenHist[index]->Fill(dlen);
	363	fChanHist[0]->Fill(type);
4c230f6d	364	if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
25c94fc8	365	if(fReadMC){
	366	if(labDs>=0) {
	367	index=GetSignalHistoIndex(iPtBin);
	368	fCosPHist[index]->Fill(cosp);
	369	fDLenHist[index]->Fill(dlen);
	370	fChanHist[1]->Fill(type);
4c230f6d	371	if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
25c94fc8	372	}else{
	373	index=GetBackgroundHistoIndex(iPtBin);
	374	fCosPHist[index]->Fill(cosp);
	375	fDLenHist[index]->Fill(dlen);
	376	fChanHist[2]->Fill(type);
4c230f6d	377	if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
25c94fc8	378	}
	379	}
	380	if(isKKpi){
	381	index=GetHistoIndex(iPtBin);
	382	Double_t invMass=d->InvMassDsKKpi();
	383	fMassHist[index]->Fill(invMass);
	384	Double_t mass01=d->InvMass2Prongs(0,1,321,321);
	385	Double_t mass12=d->InvMass2Prongs(1,2,321,211);
	386	fDalitz[index]->Fill(mass01,mass12);
4c230f6d	387	if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	388	if(fReadMC){
	389	if(labDs>=0) {
	390	index=GetSignalHistoIndex(iPtBin);
	391	fMassHist[index]->Fill(invMass);
	392	fDalitz[index]->Fill(mass01,mass12);
4c230f6d	393	if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	394	}else{
	395	index=GetBackgroundHistoIndex(iPtBin);
	396	fMassHist[index]->Fill(invMass);
	397	fDalitz[index]->Fill(mass01,mass12);
4c230f6d	398	if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	399	}
	400	}
	401	}
	402	if(ispiKK){
	403	index=GetHistoIndex(iPtBin);
	404	Double_t invMass=d->InvMassDspiKK();
	405	fMassHist[index]->Fill(invMass);
4c230f6d	406	if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	407	if(fReadMC){
	408	if(labDs>=0) {
	409	index=GetSignalHistoIndex(iPtBin);
	410	fMassHist[index]->Fill(invMass);
4c230f6d	411	if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	412	}else{
	413	index=GetBackgroundHistoIndex(iPtBin);
	414	fMassHist[index]->Fill(invMass);
4c230f6d	415	if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
25c94fc8	416	}
	417	}
	418	}
	419	}
	420	if(unsetvtx) d->UnsetOwnPrimaryVtx();
	421	}
	422
	423
	424	PostData(1,fOutput);
	425	return;
	426	}
	427
	428	//_________________________________________________________________
	429
	430	void AliAnalysisTaskSEDs::Terminate(Option_t /option*/)
	431	{
	432	// Terminate analysis
	433	//
	434	if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
	435	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	436	if (!fOutput) {
	437	printf("ERROR: fOutput not available\n");
	438	return;
	439	}
	440	fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
	441	fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
	442	fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
	443	fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
	444	fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
	445	fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
	446	fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
	447
	448
	449	TString hisname;
	450	Int_t index;
	451	for(Int_t i=0;i<fNPtBins;i++){
	452
	453	index=GetHistoIndex(i);
	454	hisname.Form("hMassAllPt%d",i);
	455	fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	456	hisname.Form("hMassAllPt%dCuts",i);
	457	fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	458	hisname.Form("hCosPAllPt%d",i);
	459	fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	460	hisname.Form("hDLenAllPt%d",i);
	461	fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	462	hisname.Form("hDalitzAllPt%d",i);
	463	fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
	464
	465	index=GetSignalHistoIndex(i);
	466	hisname.Form("hMassSigPt%d",i);
	467	fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	468	hisname.Form("hMassSigPt%dCuts",i);
	469	fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	470	hisname.Form("hCosPSigPt%d",i);
	471	fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	472	hisname.Form("hDLenSigPt%d",i);
	473	fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
	474	hisname.Form("hDalitzSigPt%d",i);
	475	fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
	476
	477	index=GetBackgroundHistoIndex(i);
	478	hisname.Form("hMassBkgPt%d",i);
	479	fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
480	hisname.Form("hMassBkgPt%dCuts",i);
481	fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
482	hisname.Form("hCosPBkgPt%d",i);
483	fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
484	hisname.Form("hDLenBkgPt%d",i);
485	fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
486	hisname.Form("hDalitzBkgPt%d",i);
487	fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
488
489	}
490	return;
491	}