#include <TDatabasePDG.h>
#include <TMath.h>
#include <TROOT.h>
+
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoDecayHF.h"
: AliAnalysisTaskSE(),
fVHFloose(0),
fVHFtight(0),
+ fReadMC(kFALSE),
fmD0PDG(),
fnbins(),
fptbins(0),
: AliAnalysisTaskSE(name),
fVHFloose(0),
fVHFtight(0),
+ fReadMC(kFALSE),
fmD0PDG(),
fnbins(),
fptbins(0),
: AliAnalysisTaskSE(name),
fVHFloose(0),
fVHFtight(0),
+ fReadMC(kFALSE),
fmD0PDG(),
fnbins(),
fptbins(0),
AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
- // In case there is an AOD handler writing a standard AOD, use the AOD
- // event in memory rather than the input (ESD) event.
- if (!aod && AODEvent() && IsStandardAOD()) aod = dynamic_cast<AliAODEvent*> (AODEvent());
+ TClonesArray *arrayD0toKpi=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();
+ arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+ }
+ } else {
+ arrayD0toKpi=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+ }
- // load D0->Kpi candidates
- TClonesArray *arrayD0toKpi =
- (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
if(!arrayD0toKpi) {
Printf("AliAnalysisTaskSECharmFraction::UserExec: D0toKpi branch not found!\n");
return;
// AOD primary vertex
AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
-
-
- // load MC particles
- TClonesArray *arrayMC =
- (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
- if(!arrayMC) {
- Printf("AliAnalysisTaskSECharmFraction::UserExec: MC particles branch not found!\n");
- return;
- }
-
- // load MC header
- AliAODMCHeader *aodmcHeader =
- (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
- if(!aodmcHeader) {
- Printf("AliAnalysisTaskSECharmFraction::UserExec: MC header branch not found!\n");
- return;
- }
-
- // MC primary vertex
+ TClonesArray *arrayMC=0;
+ AliAODMCHeader *aodmcHeader=0;
Double_t vtxTrue[3];
- aodmcHeader->GetVertex(vtxTrue);
+ if(fReadMC){
+ // load MC particles
+ arrayMC =
+ (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if(!arrayMC) {
+ Printf("AliAnalysisTaskSECharmFraction::UserExec: MC particles branch not found!\n");
+ return;
+ }
+
+ // load MC header
+ aodmcHeader =
+ (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!aodmcHeader) {
+ Printf("AliAnalysisTaskSECharmFraction::UserExec: MC header branch not found!\n");
+ return;
+ }
+
+ // MC primary vertex
+ aodmcHeader->GetVertex(vtxTrue);
+ }
if (!aod) {
Printf("ERROR: aod not available");
return;
}
+
//histogram filled with 1 for every AOD
fNentries->Fill(1);
PostData(1,fNentries);
-
+
//Printf("There are %d tracks in this event", aod->GetNumberOfTracks());
// Int_t nTotHF=0,nTotDstar=0,nTot3Prong=0;
Double_t invMassD0,invMassD0bar,ptD0,massmumtrue;
- AliAODRecoDecayHF *aodDMC=0x0;// to be used to create a fake true sec vertex
+ AliAODRecoDecayHF *aodDMC=0;// to be used to create a fake true sec vertex
// make trkIDtoEntry register (temporary)
Int_t trkIDtoEntry[100000];
for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
// cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
- if(aodDMC!=0x0)delete aodDMC;
+ if(aodDMC)delete aodDMC;
isPeakD0=kFALSE;
isPeakD0bar=kFALSE;
if((isSideBandD0||isSideBandD0bar)&&!(isPeakD0||isPeakD0bar))isSideBand=kTRUE;// TEMPORARY, NOT DONE IN THE METHOD CALLED ABOVE ONLY FOR FURTHER SIDE BAND STUDY
// INVESTIGATE SIGNAL TYPE : ACCESS TO MC INFORMATION
- aodDMC=GetD0toKPiSignalType(d,arrayMC,signallevel,massmumtrue,vtxTrue);
+ if(fReadMC){
+ aodDMC=GetD0toKPiSignalType(d,arrayMC,signallevel,massmumtrue,vtxTrue);
+ }
+ else signallevel=0;
if(!isinacceptance)signallevel=9;
fSignalType->Fill(signallevel);
// CANDIDATE VARIABLES
//NO CUTS Case: force okD0 and okD0bar = kTRUE
- if(signallevel==1)FillHistos(d,flistNoCutsSignal,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+ if(signallevel==1||signallevel==0)FillHistos(d,flistNoCutsSignal,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==2)FillHistos(d,flistNoCutsFromDstar,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==3||signallevel==4)FillHistos(d,flistNoCutsFromB,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10||signallevel==9)FillHistos(d,flistNoCutsBack,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==5||signallevel==6)FillHistos(d,flistNoCutsOther,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
//LOOSE CUTS Case
- if(signallevel==1)FillHistos(d,flistLsCutsSignal,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+ if(signallevel==1||signallevel==0)FillHistos(d,flistLsCutsSignal,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==2)FillHistos(d,flistLsCutsFromDstar,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==3||signallevel==4)FillHistos(d,flistLsCutsFromB,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10)FillHistos(d,flistLsCutsBack,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==5||signallevel==6)FillHistos(d,flistLsCutsOther,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
//TIGHT CUTS Case
- if(signallevel==1)FillHistos(d,flistTghCutsSignal,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+ if(signallevel==1||signallevel==0)FillHistos(d,flistTghCutsSignal,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==2)FillHistos(d,flistTghCutsFromDstar,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==3||signallevel==4)FillHistos(d,flistTghCutsFromB,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10)FillHistos(d,flistTghCutsBack,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
else if(signallevel==5||signallevel==6)FillHistos(d,flistTghCutsOther,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
- if(aodDMC!=0x0){
+ if(aodDMC){
delete aodDMC;
- aodDMC=0x0;
+ aodDMC=0;
}
if(unsetvtx) d->UnsetOwnPrimaryVtx();
// 11: end of the method without output
// 12: different result than MatchToMC method
- AliAODMCParticle *mum1=0x0;
- AliAODMCParticle *b1=0x0,*b2=0x0;
- AliAODMCParticle *grandmoth1=0x0;
+ AliAODMCParticle *mum1=0;
+ AliAODMCParticle *b1=0,*b2=0;
+ AliAODMCParticle *grandmoth1=0;
massMumTrue=-1;
Int_t pdgmum,dglabels[2],matchtoMC;
// get daughter AOD tracks
AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
- AliAODRecoDecayHF *aodDMC=0x0;
- if(trk0==0x0||trk1==0x0){
+ AliAODRecoDecayHF *aodDMC=0;
+ if(!trk0 || !trk1){
AliDebug(2,"Delete tracks I AM \n");
signaltype=-1;
}
dglabels[0]=trk0->GetLabel();
dglabels[1]=trk1->GetLabel();
- if(dglabels[0]==-1||dglabels[1]==-1){
+ if(dglabels[0]<0 || dglabels[1]<0){
AliDebug(2,"HERE I AM \n");
//fake tracks
b1=(AliAODMCParticle*)arrayMC->At(trk0->GetLabel());
b2=(AliAODMCParticle*)arrayMC->At(trk1->GetLabel());
- if(b1->GetMother()==-1||b2->GetMother()==-1){
+ if(!b1 || !b2) {
+ //Tracks with no mother ????? FAKE DECAY VERTEX
+
+ signaltype=10;
+ return aodDMC;
+ }
+
+ if(b1->GetMother()<0||b2->GetMother()<0){
//Tracks with no mother ????? FAKE DECAY VERTEX
-
signaltype=10;
return aodDMC;
matchtoMC=d->MatchToMC(421,arrayMC,2,pdgdaughters);
aodDMC=ConstructFakeTrueSecVtx(b1,b2,mum1,primaryVtx);
- if(aodDMC==0x0){
+ if(aodDMC){
signaltype=10;
return aodDMC;
}
}
- if(mum1->GetMother()==-1){
+ if(mum1->GetMother()<0){
// A particle coming from nothing
signaltype=10;
return aodDMC;
*/
// if(fcheckMCD0){//THIS CHECK IS NEEDED TO AVOID POSSIBLE FAILURE IN THE SECOND WHILE, FOR DEBUGGING
while(TMath::Abs(grandmoth1->GetPdgCode())!=4&&TMath::Abs(grandmoth1->GetPdgCode())!=5){
- if(grandmoth1->GetMother()==-1){
+ if(grandmoth1->GetMother()<0){
//### THE FOLLOWING IN CASE OF DEBUGGING ##########à
/*printf("mother=-1, pdgcode: %d \n",grandmoth1->GetPdgCode());
Int_t son=grandmoth1->GetDaughter(0);
AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::ConstructFakeTrueSecVtx(const AliAODMCParticle *b1, const AliAODMCParticle *b2, const AliAODMCParticle *mum,Double_t *primaryVtxTrue){
// CONSTRUCT A FAKE TRUE SECONDARY VERTEX (aodDMC)
//!!!NOTE THAT ONLY ONE MOTHER IS CONSIDERED: THE METHOD REQUIRES THE DAUGHTERS COME FROM THE SAME MOTHER !!
- if(b1==0x0||b2==0x0)return 0x0;
- if(mum==0x0)return 0x0;
+ if(!b1 || !b2)return 0;
+ if(!mum)return 0;
Double_t pD[3],xD[3],pXtrTrue[2],pYtrTrue[2],pZtrTrue[2],xtr1[3],xtr2[3];
Int_t charge[2]={0,0};
if(b1->Charge()==-1)charge[0]=1;
else {
if(b2->Charge()==-1){
//printf("Same charges for prongs \n");
- return 0x0;
+ return 0;
}
charge[1]=1;
}
pYtrTrue[charge[0]]=b1->Py();
pZtrTrue[charge[0]]=b1->Pz();
if(!b1->XvYvZv(xtr1)){
- return 0x0;
+ return 0;
}
pXtrTrue[charge[1]]=b2->Px();
if(!mum->PxPyPz(pD)){
//printf("!D from B:Get momentum failed \n");
- return 0x0;
+ return 0;
}
if(!mum->XvYvZv(xD)){
//printf("!D from B:Get position failed \n");
- return 0x0;
+ return 0;
}
/* ############ THIS HAPPENS FROM TIME TO TIME: NUMERIC PROBLEM KNOWN #################
if(pXtrTrue[0]+pXtrTrue[1]!=pD[0]){
if(!b2->XvYvZv(xtr2)){
- return 0x0;
+ return 0;
}
Double_t d0dummy[2]={0.,0.};//TEMPORARY : dummy d0 for AliAODRecoDeay constructor
AliAODRecoDecayHF* aodDMC=new AliAODRecoDecayHF(primaryVtxTrue,xD,2,0,pXtrTrue,pYtrTrue,pZtrTrue,d0dummy);
if(okD0)((TH1F*)list->FindObject(str.Data()))->Fill(invMassD0);
if(okD0bar)((TH1F*)list->FindObject(str.Data()))->Fill(invMassD0bar);
}
-
- if(massmumtrue>0.){
- str="hMassTrue";
- str.Append(namehist.Data());
- ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
-
- if(isPeakD0||isPeakD0bar){
+ if(fReadMC){
+ if(massmumtrue>0.){
str="hMassTrue";
str.Append(namehist.Data());
- str.Append("PM");
- ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
- }
- if(isSideBand){
- str="hMassTrue";
- str.Append(namehist.Data());
- str.Append("SB");
((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+
+ if(isPeakD0||isPeakD0bar){
+ str="hMassTrue";
+ str.Append(namehist.Data());
+ str.Append("PM");
+ ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+ }
+ if(isSideBand){
+ str="hMassTrue";
+ str.Append(namehist.Data());
+ str.Append("SB");
+ ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+ }
}
}
-
// ################ D0 Impact Parameter Histos #####################
if((isPeakD0&&okD0)||(isPeakD0bar&&okD0bar)){
str="hd0D0";
((TH1F*)list->FindObject(str.Data()))->Fill(d->ImpParXY()*10000.);
- if(vtxTrue){
+ if(fReadMC && vtxTrue){
str="hd0D0VtxTrue";
str.Append(namehist.Data());
str.Append("PM");
((TH1F*)list->FindObject(str.Data()))->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);
}
- if(aodDMC!=0x0){
+ if(fReadMC && aodDMC){
aodDMC->Print("");
aodDMC->ImpParXY();
aodDMC->Print("");
((TH1F*)list->FindObject(str.Data()))->Fill(d->ImpParXY()*10000.);
- if(vtxTrue){
+ if(fReadMC&&vtxTrue){
str="hd0D0VtxTrue";
str.Append(namehist.Data());
str.Append("SB");
}
- if(aodDMC!=0x0){
+ if(fReadMC && aodDMC){
str="hMCd0D0";
str.Append(namehist.Data());
str.Append("SB");
void AliAnalysisTaskSECharmFraction::SetNPtBins(Int_t nbins,const Double_t *ptbins){
- if((fptbins)!=0x0)delete fptbins;
+ // SET THE PT BINS
+ if(fptbins)delete fptbins;
fnbins=nbins;fptbins=new Double_t[fnbins];
memcpy(fptbins,ptbins,fnbins*sizeof(Double_t));
return;