--- /dev/null
+//DEFINITION OF A FEW CONSTANTS
+//
+// binning method from C.Zampolli
+//
+// general
+const Double_t ymin = -2.1 ;
+const Double_t ymax = 2.1 ;
+//soft pion
+const Double_t ptmin_0_1 = 0.0 ;
+const Double_t ptmax_0_1 = 1.0 ;
+const Double_t ptmin_1_2 = 1.0 ;
+const Double_t ptmax_1_2 = 2.0 ;
+const Double_t ptmin_2_10 = 2.0 ;
+const Double_t ptmax_2_10 = 15.0 ;
+//D0 and D0 prongs
+const Double_t ptmin_0_4 = 0.0 ;
+const Double_t ptmax_0_4 = 4.0 ;
+const Double_t ptmin_4_8 = 4.0 ;
+const Double_t ptmax_4_8 = 8.0 ;
+const Double_t ptmin_8_10 = 8.0 ;
+const Double_t ptmax_8_10 = 20.0 ;
+const Double_t cosmin = -1.05;
+const Double_t cosmax = 1.05;
+const Double_t cTmin = 0; // micron
+const Double_t cTmax = 500; // micron
+const Double_t dcamin = 0; // micron
+const Double_t dcamax = 500; // micron
+const Double_t d0min = -1000; // micron
+const Double_t d0max = 1000; // micron
+const Double_t d0xd0min = -100000; // micron
+const Double_t d0xd0max = 100000; // micron
+const Double_t phimin = 0.0;
+const Int_t mintrackrefsTPC = 2 ;
+const Int_t mintrackrefsITS = 3 ;
+const Int_t charge = 1 ;
+const Int_t minclustersTPC = 50 ;
+// cuts
+const Double_t ptmin = 0.05;
+const Double_t ptmax = 9999.;
+const Double_t etamin = -0.9;
+const Double_t etamax = 0.9;
+const Double_t zmin = -15;
+const Double_t zmax = 15;
+const Int_t minITSClusters = 3;
+const Int_t minITSClustersSoft = 2;
+//----------------------------------------------------
+
+AliCFTaskForDStarAnalysis *AddTaskCFDStar()
+{
+
+ //CONTAINER DEFINITION
+ Info("AliCFTaskForDStarAnalysis","SETUP CONTAINER");
+ //the sensitive variables, their indices
+ UInt_t ipt = 0;
+ UInt_t iy = 1;
+ UInt_t icosThetaStar = 2;
+ UInt_t ipTpi = 3;
+ UInt_t ipTD0 = 4;
+ UInt_t icT = 5;
+ UInt_t idca = 6;
+ UInt_t id0pi = 7;
+ UInt_t id0K = 8;
+ UInt_t id0xd0 = 9;
+ UInt_t ipointing = 10;
+ UInt_t iphi = 11;
+ UInt_t iz = 12;
+ UInt_t ipTD0pi = 13;
+ UInt_t ipTD0K = 14;
+
+ const Double_t phimax = 2*TMath::Pi();
+
+ //Setting up the container grid...
+ UInt_t nstep = 8; //number of selection steps
+ const Int_t nvar = 15 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi
+ const Int_t nbin0_0_4 = 8 ; //bins in pt from 0 to 4 GeV
+ const Int_t nbin0_4_8 = 4 ; //bins in pt from 4 to 8 GeV
+ const Int_t nbin0_8_10 = 2 ; //bins in pt from 8 to 10 GeV
+ const Int_t nbin1 = 30 ; //bins in y
+ const Int_t nbin2 = 30 ; //bins in cosThetaStar
+ // soft pion and D0 from D*
+ const Int_t nbin3_0_1 = 8 ; //bins in ptPi from 0 to 4 GeV
+ const Int_t nbin3_1_2 = 1 ; //bins in ptPi from 4 to 8 GeV
+ const Int_t nbin3_2_10 = 1 ; //bins in ptPi from 8 to 10 GeV
+ const Int_t nbin4_0_4 = 8 ; //bins in ptD0 from 0 to 4 GeV
+ const Int_t nbin4_4_8 = 3 ; //bins in ptD0 from 4 to 8 GeV
+ const Int_t nbin4_8_10 = 1 ; //bins in ptD0 from 8 to 10 GeV
+ // D0 prongs - cutting variables
+ const Int_t nbin5 = 20 ; //bins in cT
+ const Int_t nbin6 = 20 ; //bins in dca
+ const Int_t nbin7 = 100 ; //bins in d0pi
+ const Int_t nbin8 = 100 ; //bins in d0K
+ const Int_t nbin9 = 80 ; //bins in d0xd0
+ const Int_t nbin10 = 100 ; //bins in cosPointingAngle
+ const Int_t nbin11 = 15 ; //bins in Phi
+ const Int_t nbin12 = 60 ; //bins in z vertex
+ // D0 prongs pt and phi
+ const Int_t nbin5_0_4 = 8 ; //bins in ptPi from 0 to 4 GeV
+ const Int_t nbin5_4_8 = 4 ; //bins in ptPi from 4 to 8 GeV
+ const Int_t nbin5_8_10 = 8 ; //bins in ptPi from 8 to 10 GeV
+ const Int_t nbin6_0_4 = 8 ; //bins in ptk from 0 to 4 GeV
+ const Int_t nbin6_4_8 = 4 ; //bins in ptk from 4 to 8 GeV
+ const Int_t nbin6_8_10 = 8 ; //bins in ptk from 8 to 10 GeV
+
+ //arrays for the number of bins in each dimension
+ Int_t iBin[nvar];
+
+ iBin[0]=nbin0_0_4+nbin0_4_8+nbin0_8_10;
+ iBin[1]=nbin1;
+ iBin[2]=nbin2;
+ iBin[3]=nbin3_0_1+nbin3_1_2+nbin3_2_10;
+ iBin[4]=nbin4_0_4+nbin4_4_8+nbin4_8_10;
+ iBin[5]=nbin5;
+ iBin[6]=nbin6;
+ iBin[7]=nbin7;
+ iBin[8]=nbin8;
+ iBin[9]=nbin9;
+ iBin[10]=nbin10;
+ iBin[11]=nbin11;
+ iBin[12]=nbin12;
+ iBin[13]=nbin5_0_4+nbin5_4_8+nbin5_8_10;
+ iBin[14]=nbin6_0_4+nbin6_4_8+nbin6_8_10;
+
+ //arrays for lower bounds :
+ Double_t *binLim0 = new Double_t[iBin[0]+1];
+ Double_t *binLim1 = new Double_t[iBin[1]+1];
+ Double_t *binLim2 = new Double_t[iBin[2]+1];
+ Double_t *binLim3 = new Double_t[iBin[3]+1];
+ Double_t *binLim4 = new Double_t[iBin[4]+1];
+ Double_t *binLim5 = new Double_t[iBin[5]+1];
+ Double_t *binLim6 = new Double_t[iBin[6]+1];
+ Double_t *binLim7 = new Double_t[iBin[7]+1];
+ Double_t *binLim8 = new Double_t[iBin[8]+1];
+ Double_t *binLim9 = new Double_t[iBin[9]+1];
+ Double_t *binLim10 = new Double_t[iBin[10]+1];
+ Double_t *binLim11 = new Double_t[iBin[11]+1];
+ Double_t *binLim12 = new Double_t[iBin[12]+1];
+ Double_t *binLim13 = new Double_t[iBin[13]+1];
+ Double_t *binLim14 = new Double_t[iBin[14]+1];
+
+ // checking limits
+ if (ptmax_0_4 != ptmin_4_8) {
+ Error("AliCFTaskForDStarAnalysis","max lim 1st range != min lim 2nd range, please check!");
+ }
+ if (ptmax_4_8 != ptmin_8_10) {
+ Error("AliCFTaskForDStarAnalysis","max lim 2nd range != min lim 3rd range, please check!");
+ }
+
+ // values for bin lower bounds
+ // pt -----------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ;
+ if (binLim0[nbin0_0_4] != ptmin_4_8) {
+ Error("AliCFDStar","Calculated bin lim for pt - 1st range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ;
+ if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10) {
+ Error("AliCFDStar","Calculated bin lim for pt - 2nd range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ;
+
+ // y -----------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin + (ymax-ymin) /nbin1*(Double_t)i ;
+
+ // cosThetaStar -----------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)cosmin + (cosmax-cosmin) /nbin2*(Double_t)i ;
+
+ // Soft ptPi ---------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin3_0_1; i++) binLim3[i]=(Double_t)ptmin_0_1 + (ptmax_0_1-ptmin_0_1)/nbin3_0_1*(Double_t)i ;
+ if (binLim3[nbin3_0_1] != ptmin_1_2) {
+ Error("AliCFDStar","Calculated bin lim for ptPi - 1st range - differs from expected!");
+ }
+ for(Int_t i=0; i<=nbin3_1_2; i++) binLim3[i+nbin3_0_1]=(Double_t)ptmin_1_2 + (ptmax_1_2-ptmin_1_2)/nbin3_1_2*(Double_t)i ;
+ if (binLim3[nbin3_0_1+nbin3_1_2] != ptmin_2_10) {
+ Error("AliCFDStar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin3_2_10; i++) binLim3[i+nbin3_0_1+nbin3_1_2]=(Double_t)ptmin_2_10 + (ptmax_2_10-ptmin_2_10)/nbin3_2_10*(Double_t)i ;
+
+ // ptD0 --------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin4_0_4; i++) binLim4[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin4_0_4*(Double_t)i ;
+ if (binLim4[nbin4_0_4] != ptmin_4_8) {
+ Error("AliCFDStar","Calculated bin lim for ptKa - 1st range - differs from expected!");
+ }
+ for(Int_t i=0; i<=nbin4_4_8; i++) binLim4[i+nbin4_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin4_4_8*(Double_t)i ;
+ if (binLim4[nbin4_0_4+nbin4_4_8] != ptmin_8_10) {
+ Error("AliCFDStar","Calculated bin lim for ptKa - 2nd range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin4_8_10; i++) binLim4[i+nbin4_0_4+nbin4_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin4_8_10*(Double_t)i ;
+
+ // D0 ptPi --------------------------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin5_0_4; i++) binLim13[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin5_0_4*(Double_t)i ;
+ if (binLim13[nbin5_0_4] != ptmin_4_8) {
+ Error("AliCFDStar","Calculated bin lim for ptPi - 1st range - differs from expected!");
+ }
+ for(Int_t i=0; i<=nbin5_4_8; i++) binLim13[i+nbin5_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin5_4_8*(Double_t)i ;
+ if (binLim13[nbin5_0_4+nbin5_4_8] != ptmin_8_10) {
+ Error("AliCFDStar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin5_8_10; i++) binLim13[i+nbin5_0_4+nbin5_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin5_8_10*(Double_t)i ;
+
+ // D0 ptK ----------------------------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin6_0_4; i++) binLim14[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin6_0_4*(Double_t)i ;
+ if (binLim14[nbin6_0_4] != ptmin_4_8) {
+ Error("AliCFDStar","Calculated bin lim for ptKa - 1st range - differs from expected!");
+ }
+ for(Int_t i=0; i<=nbin6_4_8; i++) binLim14[i+nbin6_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin6_4_8*(Double_t)i ;
+ if (binLim14[nbin6_0_4+nbin6_4_8] != ptmin_8_10) {
+ Error("AliCFDStar","Calculated bin lim for ptKa - 2nd range - differs from expected!\n");
+ }
+ for(Int_t i=0; i<=nbin6_8_10; i++) binLim14[i+nbin6_0_4+nbin6_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin6_8_10*(Double_t)i ;
+
+ // cT ---------------------------------------------------------------------------------------------------------------
+ for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)cTmin + (cTmax-cTmin) /nbin5*(Double_t)i ;
+
+ // dca
+ for(Int_t i=0; i<=nbin6; i++) binLim6[i]=(Double_t)dcamin + (dcamax-dcamin) /nbin6*(Double_t)i ;
+
+ // d0pi
+ for(Int_t i=0; i<=nbin7; i++) binLim7[i]=(Double_t)d0min + (d0max-d0min) /nbin7*(Double_t)i ;
+
+ // d0K
+ for(Int_t i=0; i<=nbin8; i++) binLim8[i]=(Double_t)d0min + (d0max-d0min) /nbin8*(Double_t)i ;
+
+ // d0xd0
+ for(Int_t i=0; i<=nbin9; i++) binLim9[i]=(Double_t)d0xd0min + (d0xd0max-d0xd0min) /nbin9*(Double_t)i ;
+
+ // cosPointingAngle
+ for(Int_t i=0; i<=nbin10; i++) binLim10[i]=(Double_t)cosmin + (cosmax-cosmin) /nbin10*(Double_t)i ;
+
+ // Phi
+ for(Int_t i=0; i<=nbin11; i++) binLim11[i]=(Double_t)phimin + (phimax-phimin) /nbin11*(Double_t)i ;
+
+ // z Primary Vertex
+ for(Int_t i=0; i<=nbin12; i++) {
+ binLim12[i]=(Double_t)zmin + (zmax-zmin) /nbin12*(Double_t)i ;
+ }
+
+ //one "container" for MC
+ AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
+ //setting the bin limits
+ container -> SetBinLimits(ipt,binLim0);
+ container -> SetBinLimits(iy,binLim1);
+ container -> SetBinLimits(icosThetaStar,binLim2);
+ container -> SetBinLimits(ipTpi,binLim3);
+ container -> SetBinLimits(ipTD0,binLim4);
+ container -> SetBinLimits(icT,binLim5);
+ container -> SetBinLimits(idca,binLim6);
+ container -> SetBinLimits(id0pi,binLim7);
+ container -> SetBinLimits(id0K,binLim8);
+ container -> SetBinLimits(id0xd0,binLim9);
+ container -> SetBinLimits(ipointing,binLim10);
+ container -> SetBinLimits(iphi,binLim11);
+ container -> SetBinLimits(iz,binLim12);
+ container -> SetBinLimits(ipTD0pi,binLim13);
+ container -> SetBinLimits(ipTD0K,binLim14);
+
+ //CREATE THE CUTS -----------------------------------------------
+
+ // Gen-Level kinematic cuts
+ AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
+
+ //Particle-Level cuts:
+ AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
+ mcGenCuts->SetRequirePdgCode(413, kTRUE); // kTRUE set in order to include D*_bar
+ mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)
+
+ // Acceptance cuts:
+ AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts");
+ AliCFTrackKineCuts *kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts");
+ kineAccCuts->SetPtRange(ptmin,ptmax);
+ kineAccCuts->SetEtaRange(etamin,etamax);
+
+ // Rec-Level kinematic cuts
+ AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
+
+ AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
+
+ AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
+
+ printf("CREATE MC KINE CUTS\n");
+ TObjArray* mcList = new TObjArray(0) ;
+ mcList->AddLast(mcKineCuts);
+ mcList->AddLast(mcGenCuts);
+
+ printf("CREATE ACCEPTANCE CUTS\n");
+ TObjArray* accList = new TObjArray(0) ;
+ accList->AddLast(kineAccCuts);
+
+ printf("CREATE RECONSTRUCTION CUTS\n");
+ TObjArray* recList = new TObjArray(0) ; // not used!!
+ recList->AddLast(recKineCuts);
+ recList->AddLast(recQualityCuts);
+ recList->AddLast(recIsPrimaryCuts);
+
+ TObjArray* emptyList = new TObjArray(0);
+
+ //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
+ printf("CREATE INTERFACE AND CUTS\n");
+ AliCFManager* man = new AliCFManager() ;
+
+ man->SetParticleContainer (container);
+ man->SetParticleCutsList(0 , mcList); // MC
+ man->SetParticleCutsList(1 , accList); // Acceptance
+ man->SetParticleCutsList(2 , emptyList); // Vertex
+ man->SetParticleCutsList(3 , emptyList); // Refit
+ man->SetParticleCutsList(4 , emptyList); // AOD
+ man->SetParticleCutsList(5 , emptyList); // AOD in Acceptance
+ man->SetParticleCutsList(6 , emptyList); // AOD with required n. of ITS clusters
+ man->SetParticleCutsList(7 , emptyList); // AOD Reco cuts
+
+ // Get the pointer to the existing analysis manager via the static access method.
+ //==============================================================================
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ ::Error("AddTaskCompareHF", "No analysis manager to connect to.");
+ return NULL;
+ }
+ //CREATE THE TASK
+ printf("CREATE TASK\n");
+ // create the task
+ AliCFTaskForDStarAnalysis *task = new AliCFTaskForDStarAnalysis("AliCFTaskForDStarAnalysis");
+ task->SetMinITSClusters(minITSClusters);
+ task->SetMinITSClustersSoft(minITSClustersSoft);
+ task->SetCFManager(man); //here is set the CF manager
+
+ Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise D* cuts
+ Int_t thnDim[4];
+
+ //first half : reconstructed
+ //second half : MC
+ thnDim[0] = iBin[0];
+ thnDim[2] = iBin[0];
+ thnDim[1] = iBin[1];
+ thnDim[3] = iBin[1];
+
+ THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",4,thnDim);
+ Double_t** binEdges = new Double_t[2];
+
+ // set bin limits
+
+ binEdges[0]= binLim0;
+ binEdges[1]= binLim1;
+
+ correlation->SetBinEdges(0,binEdges[0]);
+ correlation->SetBinEdges(2,binEdges[0]);
+
+ correlation->SetBinEdges(1,binEdges[1]);
+ correlation->SetBinEdges(3,binEdges[1]);
+
+ correlation->Sumw2();
+
+ // correlation matrix ready
+ //------------------------------------------------//
+
+ task->SetCorrelationMatrix(correlation); // correlation matrix for unfolding
+
+ // Create and connect containers for input/output
+
+ // ------ input data ------
+ AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
+
+ // ----- output data -----
+
+ TString outputfile = AliAnalysisManager::GetCommonFileName();
+ outputfile += ":PWG3_D2H_CFtaskDStar";
+
+ //now comes user's output objects :
+ // output TH1I for event counting
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("CFDSchist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+ // output Correction Framework Container (for acceptance & efficiency calculations)
+ AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("CFDSccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+ // Unfolding - correlation matrix
+ AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("CFDScorr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+ mgr->AddTask(task);
+
+ mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task,1,coutput1);
+ mgr->ConnectOutput(task,2,coutput2);
+ mgr->ConnectOutput(task,3,coutput3);
+
+ return task;
+}
+
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------------
+// Class for DStar corrections:
+//
+// The D0 cutting varibles position in the container and container
+// binning method from a C.Zampolli example
+// In this way a simple comparison between D0 and D0 from D* is possible.
+//
+//-----------------------------------------------------------------------
+// Author : A.Grelli, Utrecht University
+//
+// a.grelli@uu.nl
+//-----------------------------------------------------------------------
+
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include <TH1I.h>
+#include <TStyle.h>
+#include <TFile.h>
+
+#include "AliCFTaskForDStarAnalysis.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliCFManager.h"
+#include "AliCFContainer.h"
+#include "AliLog.h"
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliESDtrack.h"
+#include "TChain.h"
+#include "THnSparse.h"
+#include "TH2D.h"
+
+//__________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis() :
+ AliAnalysisTaskSE(),
+ fCFManager(0x0),
+ fHistEventsProcessed(0x0),
+ fCorrelation(0x0),
+ fCountRecoDStarSel(0),
+ fEvents(0),
+ fMinITSClusters(5),
+ fMinITSClustersSoft(4),
+ fAcceptanceUnf(kTRUE)
+{
+ //
+ //Default ctor
+ //
+}
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const Char_t* name) :
+ AliAnalysisTaskSE(name),
+ fCFManager(0x0),
+ fHistEventsProcessed(0x0),
+ fCorrelation(0x0),
+ fCountRecoDStarSel(0),
+ fEvents(0),
+ fMinITSClusters(5),
+ fMinITSClustersSoft(4),
+ fAcceptanceUnf(kTRUE)
+{
+ //
+ // Constructor. Initialization of Inputs and Outputs
+ //
+ Info("AliCFTaskForDStarAnalysis","Calling Constructor");
+
+ DefineOutput(1,TH1I::Class());
+ DefineOutput(2,AliCFContainer::Class());
+ DefineOutput(3,THnSparseD::Class());
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis& AliCFTaskForDStarAnalysis::operator=(const AliCFTaskForDStarAnalysis& c)
+{
+ //
+ // Assignment operator
+ //
+ if (this!=&c) {
+ AliAnalysisTaskSE::operator=(c) ;
+ fCFManager = c.fCFManager;
+ fHistEventsProcessed = c.fHistEventsProcessed;
+ }
+ return *this;
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const AliCFTaskForDStarAnalysis& c) :
+ AliAnalysisTaskSE(c),
+ fCFManager(c.fCFManager),
+ fHistEventsProcessed(c.fHistEventsProcessed),
+ fCorrelation(c.fCorrelation),
+ fCountRecoDStarSel(c.fCountRecoDStarSel),
+ fEvents(c.fEvents),
+ fMinITSClusters(c.fMinITSClusters),
+ fMinITSClustersSoft(c.fMinITSClustersSoft),
+ fAcceptanceUnf(c.fAcceptanceUnf)
+{
+ //
+ // Copy Constructor
+ //
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
+ //
+ //destructor
+ //
+ if (fCFManager) delete fCFManager ;
+ if (fHistEventsProcessed) delete fHistEventsProcessed ;
+ if (fCorrelation) delete fCorrelation ;
+}
+
+//_________________________________________________
+void AliCFTaskForDStarAnalysis::UserExec(Option_t *)
+{
+ //
+ // Main loop function
+ //
+
+ if (!fInputEvent) {
+ Error("UserExec","NO EVENT FOUND!");
+ return;
+ }
+
+ fEvents++;
+
+ if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
+ 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("Dstar");
+ }
+ } else {
+ arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
+ }
+
+ if (!arrayDStartoD0pi) {
+ AliError("Could not find array of HF vertices");
+ return;
+ }
+
+ fCFManager->SetRecEventInfo(aodEvent);
+ fCFManager->SetMCEventInfo(aodEvent);
+
+ // event selection
+ Double_t containerInput[14] ;
+ Double_t containerInputMC[14] ;
+
+ //loop on the MC event
+
+ TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) {
+ AliError("Could not find Monte-Carlo in AOD");
+ return;
+ }
+
+ AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+ if (!mcHeader) {
+ AliError("Could not find MC Header in AOD");
+ return;
+ }
+
+ // AOD primary vertex
+ AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+ Double_t zPrimVertex = vtx1->GetZ();
+ Double_t zMCVertex = mcHeader->GetVtxZ();
+ Bool_t vtxFlag = kTRUE;
+ TString title=vtx1->GetTitle();
+ if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
+
+ for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
+ AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+ if (!mcPart) {
+ AliWarning("MC Particle not found in tree, skipping");
+ continue;
+ }
+
+ // check the MC-level cuts
+ if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
+
+ // fill the container for Gen-level selection
+ Double_t vectorMC[9] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
+
+ if (GetDStarMCParticle(mcPart, mcArray, vectorMC)){
+
+ containerInputMC[0] = vectorMC[0];
+ containerInputMC[1] = vectorMC[1] ;
+ containerInputMC[2] = vectorMC[2] ;
+ containerInputMC[3] = vectorMC[3] ;
+ containerInputMC[4] = vectorMC[4] ;
+ containerInputMC[5] = vectorMC[5] ;
+ containerInputMC[6] = 0.;
+ containerInputMC[7] = 0.;
+ containerInputMC[8] = 0.;
+ containerInputMC[9] = -100000.;
+ containerInputMC[10] = 1.01;
+ containerInputMC[11] = vectorMC[6];
+ containerInputMC[12] = zMCVertex; // z vertex
+ containerInputMC[13] = vectorMC[7]; // pt D0 pion
+ containerInputMC[14] = vectorMC[8]; // pt D0 kaon
+
+ fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
+
+ // check the MC-Acceptance level cuts
+ // since standard CF functions are not applicable, using Kine Cuts on daughters
+
+ Int_t daughter0 = mcPart->GetDaughter(0);
+ Int_t daughter1 = mcPart->GetDaughter(1);
+
+ AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
+
+ //D0
+ AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
+ // Soft Pion
+ AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
+
+ // Acceptance variables for the soft pion
+ Double_t eta1 = mcPartDaughter1->Eta();
+ Double_t pt1 = mcPartDaughter1->Pt();
+
+ Int_t daughD00 = 0;
+ Int_t daughD01 = 0;
+
+ // Just to be sure to take the right particles
+ if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
+ daughD00 = mcPartDaughter0->GetDaughter(0);
+ daughD01 = mcPartDaughter0->GetDaughter(1);
+ }else{
+ daughD00 = mcPartDaughter1->GetDaughter(0);
+ daughD01 = mcPartDaughter1->GetDaughter(1);
+ }
+
+ // the D0 daughters
+ AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
+ AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
+
+ if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
+ AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
+ }
+
+ // D0 daughters - needed for acceptance
+ Double_t theD0pt0 = mcD0PartDaughter0->Pt();
+ Double_t theD0pt1 = mcD0PartDaughter1->Pt();
+ Double_t theD0eta0 = mcD0PartDaughter0->Eta();
+ Double_t theD0eta1 = mcD0PartDaughter1->Eta();
+
+ // ACCEPTANCE REQUESTS ---------
+
+ // soft pion
+ Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
+ // Do daughters
+ Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1);
+ Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
+
+ if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
+
+ AliDebug(2, "D* Daughter particles in acceptance");
+
+ fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
+
+ // check on the vertex
+ if (vtxFlag){
+ // filling the container if the vertex is ok
+ fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
+
+ Bool_t refitFlag = kTRUE;
+ for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
+ AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
+
+ // refit only for D0 daughters
+ if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
+ if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
+ refitFlag = kFALSE;
+ }
+ }
+ }
+ if (refitFlag){ // refit
+
+ fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
+
+ } // end of refit
+ } // end of vertex
+ } //end of acceptance
+ } // end of MC D*
+ else {
+ AliDebug(3,"Problems in filling the container");
+ continue;
+ }
+ } // end of MC loop
+
+ //rec
+ AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
+
+ //D* and D0 prongs needed to MatchToMC method
+ Int_t pdgDgDStartoD0pi[2]={421,211};
+ Int_t pdgDgD0toKpi[2]={321,211};
+
+ for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
+
+ // D* candidates
+ AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
+
+ // D0 from the reco cascade
+ AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
+ Bool_t unsetvtx=kFALSE;
+
+ // needed for pointing angle
+ if(!D0particle->GetOwnPrimaryVtx()) {
+ D0particle->SetOwnPrimaryVtx(vtx1);
+ unsetvtx=kTRUE;
+ }
+
+ // find associated MC particle for D* ->D0toKpi
+ Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
+
+ // find D0->Kpi ... needed in the following
+ Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
+
+ if (mcLabel == -1 || mcD0Label ==-1)
+ {
+ AliDebug(2,"No MC particle found");
+ }
+ else {
+
+ // the D* and the D0 in MC
+ AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
+ AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
+
+ if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
+ AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
+ continue;
+ }
+
+ // soft pion
+ AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
+
+ //D0tokpi
+ AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
+ AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
+
+ // check if associated MC v0 passes the cuts
+ if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) {
+ AliDebug(2, "Skipping the particles due to cuts");
+ continue;
+ }
+
+ // fill the container...
+ Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
+ Double_t rapidity = dstarD0pi->YDstar();
+ Double_t cosThetaStar = 9999.;
+ Double_t pTpi = 0.;
+ Double_t pTK = 0.;
+ Double_t dca = (D0particle->GetDCA())*1E4;
+ Double_t d0pi = 0.;
+ Double_t d0K = 0.;
+ Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
+ Double_t cosPointingAngle = D0particle->CosPointingAngle();
+ Double_t phi = D0particle->Phi();
+
+ // Select D0 cutting variables
+ Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
+
+ // D0 related variables
+ if (pdgCode > 0){
+
+ cosThetaStar = D0particle->CosThetaStarD0();
+ pTpi = D0particle->PtProng(0);
+ pTK = D0particle->PtProng(1);
+ d0pi = (D0particle->Getd0Prong(0))*1E4;
+ d0K = (D0particle->Getd0Prong(1))*1E4;
+
+ }
+ else {
+
+ cosThetaStar = D0particle->CosThetaStarD0bar();
+ pTpi = D0particle->PtProng(1);
+ pTK = D0particle->PtProng(0);
+ d0pi = (D0particle->Getd0Prong(1))*1E4;
+ d0K = (D0particle->Getd0Prong(0))*1E4;
+
+ }
+
+ // ct of the D0 from D*
+ Double_t cT = D0particle->CtD0();
+
+ containerInput[0] = pt;
+ containerInput[1] = rapidity;
+ containerInput[2] = cosThetaStar;
+ containerInput[3] = track2->Pt();
+ containerInput[4] = D0particle->Pt();
+ containerInput[5] = cT*1.E4; // in micron
+ containerInput[6] = dca; // in micron
+ containerInput[7] = d0pi; // in micron
+ containerInput[8] = d0K; // in micron
+ containerInput[9] = d0xd0; // in micron^2
+ containerInput[10] = cosPointingAngle; // in micron
+ containerInput[11] = phi;
+ containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
+ containerInput[13] = pTpi; // D0 pion
+ containerInput[14] = pTK; // D0 kaon
+
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
+
+ // refit in ITS and TPC for D0 daughters
+ if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
+ continue;
+ }
+
+ // reft in ITS for soft pion
+ if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
+ continue;
+ }
+
+ // cut in acceptance for the soft pion and for the D0 daughters
+ Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
+ Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
+
+ // soft pion acceptance ... is it fine 0.9?????
+ Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
+
+ if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
+ AliDebug(2,"D* reco daughters in acceptance");
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
+
+ if(fAcceptanceUnf){
+ Double_t fill[4]; //fill response matrix
+
+ // dimensions 0&1 : pt,eta (Rec)
+ fill[0] = pt ;
+ fill[1] = rapidity;
+ // dimensions 2&3 : pt,eta (MC)
+ fill[2] = mcVtxHFDStar->Pt();
+ fill[3] = mcVtxHFDStar->Y();
+ fCorrelation->Fill(fill);
+ }
+
+ // cut on the min n. of clusters in ITS for the D0 and soft pion
+ Int_t ncls0=0,ncls1=0,ncls2=0;
+ for(Int_t l=0;l<6;l++) {
+ if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
+ if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
+ if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
+ }
+ // see AddTask for soft pion ITS clusters request
+ AliDebug(2, Form("n clusters = %d", ncls0));
+
+ if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
+
+ // D0 cuts optimized for D* analysis
+ Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
+
+ // needed for cuts
+ Double_t theD0pt = D0particle->Pt();
+
+ if (theD0pt <= 1){ // first bin not optimized
+ cuts[0] = 400;
+ cuts[1] = 0.8;
+ cuts[2] = 0.21;
+ cuts[3] = 500;
+ cuts[4] = 500;
+ cuts[5] = -20000;
+ cuts[6] = 0.6;
+ }
+ else if (theD0pt > 1 && theD0pt <= 2){
+ cuts[0] = 200;
+ cuts[1] = 0.7;
+ cuts[2] = 0.8;
+ cuts[3] = 210;
+ cuts[4] = 210;
+ cuts[5] = -20000;
+ cuts[6] = 0.9;
+ }
+ else if (theD0pt > 2 && theD0pt <= 3){
+ cuts[0] = 400;
+ cuts[1] = 0.8;
+ cuts[2] = 0.8;
+ cuts[3] = 420;
+ cuts[4] = 350;
+ cuts[5] = -8500;
+ cuts[6] = 0.9;
+ }
+ else if (theD0pt > 3 && theD0pt <= 5){
+ cuts[0] = 160;
+ cuts[1] = 1.0;
+ cuts[2] = 1.2;
+ cuts[3] = 560;
+ cuts[4] = 420;
+ cuts[5] = -8500;
+ cuts[6] = 0.9;
+ }
+ else if (theD0pt > 5){
+ cuts[0] = 800;
+ cuts[1] = 1.0;
+ cuts[2] = 1.2;
+ cuts[3] = 700;
+ cuts[4] = 700;
+ cuts[5] = 10000;
+ cuts[6] = 0.9;
+ }
+ if (dca < cuts[0]
+ && TMath::Abs(cosThetaStar) < cuts[1]
+ && pTpi > cuts[2]
+ && pTK > cuts[2]
+ && TMath::Abs(d0pi) < cuts[3]
+ && TMath::Abs(d0K) < cuts[4]
+ && d0xd0 < cuts[5]
+ && cosPointingAngle > cuts[6]
+ ){
+
+ AliDebug(2,"Particle passed D* selection cuts");
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;
+
+ if(!fAcceptanceUnf){ // unfolding
+
+ Double_t fill[4]; //fill response matrix
+
+ // dimensions 0&1 : pt,eta (Rec)
+ fill[0] = pt ;
+ fill[1] = rapidity;
+ // dimensions 2&3 : pt,eta (MC)
+ fill[2] = mcVtxHFDStar->Pt();
+ fill[3] = mcVtxHFDStar->Y();
+
+ fCorrelation->Fill(fill);
+ }
+ }
+ }
+ }
+ }
+ if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
+ } // end loop on D*->Kpipi
+
+ fHistEventsProcessed->Fill(0);
+
+ PostData(1,fHistEventsProcessed) ;
+ PostData(2,fCFManager->GetParticleContainer()) ;
+ PostData(3,fCorrelation) ;
+}
+
+
+//___________________________________________________________________________
+void AliCFTaskForDStarAnalysis::Terminate(Option_t*)
+{
+ // The Terminate()
+ AliAnalysisTaskSE::Terminate();
+
+ // draw correlation matrix
+ AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+ if(!cont) {
+ printf("CONTAINER NOT FOUND\n");
+ return;
+ }
+}
+
+//___________________________________________________________________________
+void AliCFTaskForDStarAnalysis::UserCreateOutputObjects() {
+
+ //
+ // useroutput
+ //
+
+ Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
+
+ //slot #1
+ OpenFile(1);
+ fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
+}
+
+//________________________________________________________________________________
+Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
+
+ //
+ // fill the D* and D0 MC container
+ //
+
+ Bool_t isDStar = kFALSE;
+
+ if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
+
+ // getting the daughters
+ Int_t daughter0 = mcPart->GetDaughter(0);
+ Int_t daughter1 = mcPart->GetDaughter(1);
+
+ AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
+ if (daughter0 == 0 || daughter1 == 0) {
+ AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
+ return isDStar;
+ }
+
+ if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
+ AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
+ return isDStar;
+ }
+
+ AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
+ AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
+ if (!mcPartDaughter0 || !mcPartDaughter1) {
+ AliWarning("D*: At least one Daughter Particle not found in tree, skipping");
+ return isDStar;
+ }
+
+ if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
+ TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
+ !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
+ TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
+ AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
+ return isDStar;
+ }
+
+ Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
+ Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
+
+ // getting vertex from daughters
+ mcPartDaughter0->XvYvZv(vtx2daughter0);
+ mcPartDaughter1->XvYvZv(vtx2daughter1);
+
+ // check if the secondary vertex is the same for both
+ if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
+ AliError("The D* daughters have different secondary vertex, skipping the track");
+ return isDStar;
+ }
+
+ AliAODMCParticle* neutralDaugh = mcPartDaughter0;
+
+ Double_t VectorD0[2] ={0.,0.};
+
+ if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
+ AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+ return isDStar;
+ }
+ // get the pT of the daughters
+
+ Double_t pTpi = 0.;
+ Double_t pTD0 = 0.;
+
+ if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
+ pTpi = mcPartDaughter0->Pt();
+ pTD0 = mcPartDaughter1->Pt();
+ }
+ else {
+ pTpi = mcPartDaughter1->Pt();
+ pTD0 = mcPartDaughter0->Pt();
+ }
+
+ vectorMC[0] = mcPart->Pt();
+ vectorMC[1] = mcPart->Y() ;
+ vectorMC[2] = 0;
+ vectorMC[3] = pTpi ;
+ vectorMC[4] = pTD0 ;
+ vectorMC[5] = 0;
+ vectorMC[6] = mcPart->Phi() ;
+ vectorMC[7] = VectorD0[0] ;
+ vectorMC[8] = VectorD0[1] ;
+
+ isDStar = kTRUE;
+
+ return isDStar;
+}
+//________________________________________________________________________________________________
+
+Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
+
+ //
+ // chack wether D0 is decaing into kpi
+ //
+
+ Bool_t isHadronic = kFALSE;
+
+ Int_t daughterD00 = neutralDaugh->GetDaughter(0);
+ Int_t daughterD01 = neutralDaugh->GetDaughter(1);
+
+ AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
+ if (daughterD00 == 0 || daughterD01 == 0) {
+ AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+ return isHadronic;
+ }
+
+ if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
+ AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
+ return isHadronic;
+ }
+
+ AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
+ AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
+ if (!mcPartDaughterD00 || !mcPartDaughterD01) {
+ AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
+ return isHadronic;
+ }
+
+ if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
+ TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
+ !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
+ TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
+ AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
+ return isHadronic;
+ }
+
+ Double_t pTD0pi = 0;
+ Double_t pTD0K = 0;
+
+
+ if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
+ pTD0pi = mcPartDaughterD00->Pt();
+ pTD0K = mcPartDaughterD01->Pt();
+ }
+ else {
+ pTD0pi = mcPartDaughterD01->Pt();
+ pTD0K = mcPartDaughterD00->Pt();
+ }
+
+ isHadronic = kTRUE;
+
+ VectorD0[0] = pTD0pi;
+ VectorD0[1] = pTD0K;
+
+ return isHadronic;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff