[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 "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSEDs.h"

ClassImp(AliAnalysisTaskSEDs)


//________________________________________________________________________
AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
  AliAnalysisTaskSE(),
  fOutput(0), 
  fHistNEvents(0),
  fReadMC(kFALSE),
  fNPtBins(0),
  fMassRange(0.2),
  fVHF(0)
{
  // Default constructor
}

//________________________________________________________________________
AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name):
  AliAnalysisTaskSE(name),
  fOutput(0),
  fHistNEvents(0),
  fReadMC(kFALSE),
  fNPtBins(0),
  fMassRange(0.2),
  fVHF(0)
{
  // Default constructor
  // Output slot #1 writes into a TList container
  Double_t ptlim[6]={0.,1.,3.,5.,10.,9999999.};  
  SetPtBins(5,ptlim);
 
  DefineOutput(1,TList::Class());  //My private output

}

//________________________________________________________________________
void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Double_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 (fVHF) {
    delete fVHF;
    fVHF = 0;
  }
}  

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

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

  gROOT->LoadMacro("ConfigVertexingHF.C");

  fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
  fVHF->PrintStatus();

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