#include <TList.h>
#include <TH1F.h>
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODVertex.h"
+#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoDecayHF4Prong.h"
#include "AliAODRecoCascadeHF.h"
#include "AliAnalysisVertexingHF.h"
#include "AliAnalysisTaskSE.h"
fOutput(0),
fNtupleCmp(0),
fHistMass(0),
+fHistNEvents(0),
fVHF(0)
{
// Default constructor
fOutput(0),
fNtupleCmp(0),
fHistMass(0),
+fHistNEvents(0),
fVHF(0)
{
// Standard constructor
fHistMass->SetMinimum(0);
fOutput->Add(fHistMass);
- fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
+ fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
+ fHistNEvents->Sumw2();
+ fHistNEvents->SetMinimum(0);
+ fOutput->Add(fHistNEvents);
+
+ OpenFile(2); // 2 is the slot number of the ntuple
+ fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
return;
}
AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
-
- // load HF vertices
- TClonesArray *inputArrayVertices =
- (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
+ TClonesArray *inputArrayVertices = 0;
+ TClonesArray *inputArrayD0toKpi = 0;
+ TClonesArray *inputArrayDstar = 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();
+ // load HF vertices
+ inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
+ // load D0->Kpi candidates
+ inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+ // load D*+ candidates
+ inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
+ }
+ } else {
+ // load HF vertices
+ inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
+ // load D0->Kpi candidates
+ inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+ // load D*+ candidates
+ inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
+ }
+
+
if(!inputArrayVertices) {
printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
return;
}
-
- // load D0->Kpi candidates
- TClonesArray *inputArrayD0toKpi =
- (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
if(!inputArrayD0toKpi) {
printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
return;
}
-
-
- // load D*+ candidates
- TClonesArray *inputArrayDstar =
- (TClonesArray*)aod->GetList()->FindObject("Dstar");
if(!inputArrayDstar) {
printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
return;
}
+ fHistNEvents->Fill(0); // count event
+ // Post the data already here
+ PostData(1,fOutput);
+
// AOD primary vertex
AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
//vtx1->Print();
Int_t nprongs,lab,okD0,okD0bar,pdg;
Bool_t unsetvtx;
- Double_t invmass,cov[6],errx,erry,errz;
+ Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
AliAODRecoDecayHF2Prong *d2=0;
AliAODRecoDecayHF3Prong *d3=0;
AliAODRecoDecayHF4Prong *d4=0;
+ Int_t pdgDgD0toKpi[2]={321,211};
+ Int_t pdgDgDplustoKpipi[3]={321,211,211};
+ Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
+
// loop over vertices
Int_t nVertices = inputArrayVertices->GetEntriesFast();
- printf("Number of vertices: %d\n",nVertices);
+ if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
+ vtx->GetXYZ(posRec);
+ vtx->GetCovarianceMatrix(covRec);
+ errx=1.; erry=1.; errz=1.;
+ if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
+ if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
+ if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
+
+
// get parent AliAODRecoDecayHF
nprongs= vtx->GetNDaughters();
+ // check that the daughters have kITSrefit and kTPCrefit
+ Bool_t allDgOK=kTRUE;
+ for(Int_t idg=0; idg<nprongs; idg++) {
+ AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
+ if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
+ if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
+ }
+ if(!allDgOK) continue;
+
switch(nprongs) {
case 2: // look for D0->Kpi
d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
- lab = d2->MatchToMC(421,mcArray);
+ if(d2->IsLikeSign()) continue;
+ if(d2->Charge() != 0) continue; // these are D*
+ lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
if(lab>=0) {
unsetvtx=kFALSE;
if(!d2->GetOwnPrimaryVtx()) {
unsetvtx=kTRUE;
}
okD0=0; okD0bar=0;
- if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
- AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
- pdg = part->GetPdgCode();
+ if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // beware: cuts may bias the resolution!
+ AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
+ pdg = dMC->GetPdgCode();
invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
+ // get a daughter for true pos of decay vertex
+ AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
+ dg0MC->XvYvZv(posTrue);
fHistMass->Fill(invmass);
// Post the data already here
PostData(1,fOutput);
-
- AliAODVertex *secVtx = d2->GetSecondaryVtx();
- secVtx->GetCovarianceMatrix(cov);
- errx=1.; erry=1.; errz=1.;
- if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
- if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
- if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
-
- fNtupleCmp->Fill(pdg,nprongs,d2->Xv(),part->Xv(),errx,d2->Yv(),part->Yv(),erry,d2->Zv(),part->Zv(),errz,d2->GetReducedChi2(),d2->Pt(),invmass);
+ Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+ (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+ (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+ (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+ (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
+ (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
+ fNtupleCmp->Fill(tmp);
PostData(2,fNtupleCmp);
}
if(unsetvtx) d2->UnsetOwnPrimaryVtx();
break;
case 3: // look for D+
d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
- lab = d3->MatchToMC(411,mcArray);
+ if(d3->IsLikeSign()) continue;
+ lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
if(lab>=0) {
unsetvtx=kFALSE;
if(!d3->GetOwnPrimaryVtx()) {
d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
unsetvtx=kTRUE;
}
- if(d3->SelectDplus(fVHF->GetDplusCuts())) {
- AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
- pdg = part->GetPdgCode();
- invmass = d3->InvMassDplus();
- AliAODVertex *secVtx = d3->GetSecondaryVtx();
- secVtx->GetCovarianceMatrix(cov);
- errx=1.; erry=1.; errz=1.;
- if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
- if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
- if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
-
- fNtupleCmp->Fill(pdg,nprongs,d3->Xv(),part->Xv(),errx,d3->Yv(),part->Yv(),erry,d3->Zv(),part->Zv(),errz,d3->GetReducedChi2(),d3->Pt(),invmass);
- PostData(2,fNtupleCmp);
- }
+ //if(d3->SelectDplus(fVHF->GetDplusCuts())) {
+ AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
+ pdg = dMC->GetPdgCode();
+ invmass = d3->InvMassDplus();
+ // get a daughter for true pos of decay vertex
+ AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
+ dg0MC->XvYvZv(posTrue);
+ Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+ (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+ (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+ (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+ (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
+ (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
+ fNtupleCmp->Fill(tmp);
+ PostData(2,fNtupleCmp);
+ //}
if(unsetvtx) d3->UnsetOwnPrimaryVtx();
}
break;
case 4: // look for D0->Kpipipi
d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
- lab = d4->MatchToMC(421,mcArray);
+ if(d4->IsLikeSign()) continue;
+ lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
if(lab>=0) {
unsetvtx=kFALSE;
if(!d4->GetOwnPrimaryVtx()) {
unsetvtx=kTRUE;
}
okD0=0; okD0bar=0;
- if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
- AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
- pdg = part->GetPdgCode();
- //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
- invmass = 10.;
-
- AliAODVertex *secVtx = d4->GetSecondaryVtx();
- secVtx->GetCovarianceMatrix(cov);
- errx=1.; erry=1.; errz=1.;
- if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
- if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
- if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
-
- fNtupleCmp->Fill(pdg,nprongs,d4->Xv(),part->Xv(),errx,d4->Yv(),part->Yv(),erry,d4->Zv(),part->Zv(),errz,d4->GetReducedChi2(),d4->Pt(),invmass);
- PostData(2,fNtupleCmp);
-
- }
+ //if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
+ AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
+ pdg = dMC->GetPdgCode();
+ //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
+ invmass = 10.; // dummy
+ // get a daughter for true pos of decay vertex
+ AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
+ dg0MC->XvYvZv(posTrue);
+ Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+ (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+ (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+ (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+ (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
+ (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
+ fNtupleCmp->Fill(tmp);
+ PostData(2,fNtupleCmp);
+ //}
if(unsetvtx) d4->UnsetOwnPrimaryVtx();
}
break;
}
fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
+ fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
- fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
+ //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
return;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff