#include "AliAODEvent.h"
+#include "AliAODRecoDecay.h"
#include "AliAODTrack.h"
#include "AliLog.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliMCEvent.h"
+#include "AliAODMCParticle.h"
+
+
#include "AliAnalysisTaskAODFilterBitQA.h"
//________________________________________________________________________
AliAnalysisTaskAODFilterBitQA::AliAnalysisTaskAODFilterBitQA(const char *name)
: AliAnalysisTaskSE(name),
+ fArrayMC(0x0),
+ fListQA(0x0),
+ useCentrality(kFALSE),
+ fillOnlySecondaries(kFALSE),
+ fillHFVertexingTracks(kFALSE),
+ fHFBranchName("D0toKpi"),
+ fBitIgnore1(-1),
+ fBitIgnore2(-1),
+ fCentralityPercentileMin(0.),
+ fCentralityPercentileMax(80.),
+ fPtMin(0),
+ fPtMax(1000),
+ fEtaMin(-10),
+ fEtaMax(10),
fHistTrackStats(0)
{
- for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
- fHistKinematics[iTrackBit] = NULL;
- fHistDCA[iTrackBit] = NULL;
- fHistDCAprop[iTrackBit] = NULL;
- fHistChiClus[iTrackBit] = NULL;
+ for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
+ for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
+ fHistKinematics[iCharge][iTrackBit] = NULL;
+ fHistDCAconstrained[iCharge][iTrackBit] = NULL;
+ fHistDCAglobal[iCharge][iTrackBit] = NULL;
+ fHistChiClus[iCharge][iTrackBit] = NULL;
+ }
}
DefineInput(0, TChain::Class());
fHistTrackStats = new TH2D("fHistTrackStats","Track statistics;Centrality;TrackFilterBit;N_{events}",100,0,100,gBitMax,0,gBitMax);
fListQA->Add(fHistTrackStats);
- for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
- fHistKinematics[iTrackBit] = new TH3D(Form("Bit%d_Kinematics",iTrackBit),Form("Bit%d_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
- fHistDCA[iTrackBit] = new TH2D(Form("Bit%d_DCA",iTrackBit),Form("Bit%d_DCA;DCA XY [AODtrack] (cm);DCA Z [AODtrack] (cm)",iTrackBit),100,-5.0,5.0,220,-11.0,11.0);
- fHistDCAprop[iTrackBit] = new TH2D(Form("Bit%d_DCAprop",iTrackBit),Form("Bit%d_DCAprop;DCA XY [Propagated] (cm);DCA Z [Propagated] (cm)",iTrackBit),100,-5.0,5.0,220,-11.0,11.0);
- fHistChiClus[iTrackBit] = new TH2D(Form("Bit%d_ChiClus",iTrackBit),Form("Bit%d_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit),100,-1.0,5.0,160,0,160.0);
- fListQA->Add(fHistKinematics[iTrackBit]);
- fListQA->Add(fHistDCA[iTrackBit]);
- fListQA->Add(fHistDCAprop[iTrackBit]);
- fListQA->Add(fHistChiClus[iTrackBit]);
+ TString sCharge[gNCharge] = {"Plus","Minus"};
+
+ for(Int_t iCharge = 0; iCharge < gNCharge; iCharge++){
+ for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
+ fHistKinematics[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_Kinematics",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_Kinematics;#eta;#varphi (rad);p_{T} (GeV/c)",iTrackBit,sCharge[iCharge].Data()),100,-1.0,1.0,100,0,TMath::Pi()*2,100,0,10);
+ fHistDCAconstrained[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_DCAconstrained",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAconstrained;DCA XY [Constrained] (cm);DCA Z [Constrained] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0);
+ fHistDCAglobal[iCharge][iTrackBit] = new TH3D(Form("Bit%d_%s_DCAglobal",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_DCAglobal;DCA X [Global] (cm);DCA Y [Global] (cm);DCA Z [Global] (cm)",iTrackBit,sCharge[iCharge].Data()),100,-5.0,5.0,100,-5.0,5.0,100,-5.0,5.0);
+fHistChiClus[iCharge][iTrackBit] = new TH2D(Form("Bit%d_%s_ChiClus",iTrackBit,sCharge[iCharge].Data()),Form("Bit%d_%s_ChiClus;#chi^{2} [Fit];N_{clus} [TPC]",iTrackBit,sCharge[iCharge].Data()),100,-1.0,5.0,160,0,160.0);
+ fListQA->Add(fHistKinematics[iCharge][iTrackBit]);
+ fListQA->Add(fHistDCAconstrained[iCharge][iTrackBit]);
+ fListQA->Add(fHistDCAglobal[iCharge][iTrackBit]);
+ fListQA->Add(fHistChiClus[iCharge][iTrackBit]);
+ }
}
// Post output data.
return;
}
-
-
+ // MC information (set if available)
+ fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
+
// check event cuts
Double_t lMultiplicityVar = -1;
if((lMultiplicityVar = IsEventAccepted(event)) < 0){
return;
}
-
- // get the accepted tracks in main event
- GetAcceptedTracks(event,lMultiplicityVar);
-
+ // fill HF vertexing tracks
+ if(fillHFVertexingTracks){
+ GetAcceptedHFVertexingTracks(event,lMultiplicityVar);
+ }
+ else{
+ // get the accepted tracks in main event
+ GetAcceptedTracks(event,lMultiplicityVar);
+ }
}
//________________________________________________________________________
// Checks the Event cuts
// still hard coded
- Double_t fCentralityPercentileMin = 0.;
- Double_t fCentralityPercentileMax = 80.;
Double_t fVxMax = 0.5;
Double_t fVyMax = 0.5;
Double_t fVzMax = 10.0;
if(TMath::Abs(vertex->GetX()) < fVxMax) {
if(TMath::Abs(vertex->GetY()) < fVyMax) {
if(TMath::Abs(vertex->GetZ()) < fVzMax) {
+
+ // get the reference multiplicty or centrality (if required)
+ if(useCentrality){
- // get the reference multiplicty or centrality
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();
- gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
-
- if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();
+ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
- return gCentrality;
-
- }//centrality range
+ if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
+ return gCentrality;
+ }
+ else{
+ return -1;
+ }//centrality range
+ }//use centrality
+
+ // if not using centrality/multiplicty, return 1
+ return 1;
+
}//Vz cut
}//Vy cut
}//Vx cut
// Fills QA histograms
- Short_t vCharge;
Double_t vEta;
- Double_t vY;
Double_t vPhi;
Double_t vPt;
- Double_t vDCAxy;
- Double_t vDCAz;
- Double_t vDCApropxy;
- Double_t vDCApropz;
+ Double_t vDCAconstrainedxy;
+ Double_t vDCAconstrainedz;
+ Double_t vDCAglobalx;
+ Double_t vDCAglobaly;
+ Double_t vDCAglobalz;
Double_t vChi2;
Double_t vClus;
- //for propagation to DCA
+ Double_t pos[3];
+ Double_t v[3];
+
const AliVVertex *vertex = event->GetPrimaryVertex();
- Double_t field = event->GetMagneticField();
- Double_t b[2];
- Double_t bCov[3];
+ vertex->GetXYZ(v);
// Loop over tracks in event
for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
continue;
}
+ // get MC information (if available)
+ if(fArrayMC && fillOnlySecondaries){
+
+ Int_t label = aodTrack->GetLabel();
+ AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));
+
+ if(mcTrack->IsPhysicalPrimary())
+ continue;
+ }
+
// track parameters
- vCharge = aodTrack->Charge();
vEta = aodTrack->Eta();
- vY = aodTrack->Y();
vPhi = aodTrack->Phi();// * TMath::RadToDeg();
vPt = aodTrack->Pt();
- vDCAxy = aodTrack->DCA();
- vDCAz = aodTrack->ZAtDCA();
+ vDCAconstrainedxy = aodTrack->DCA();
+ vDCAconstrainedz = aodTrack->ZAtDCA();
vChi2 = aodTrack->Chi2perNDF();
vClus = aodTrack->GetTPCNcls();
- // propagate again to DCA, use copy to propagate (in order not to overwrite track parameters)
- AliAODTrack copy(*aodTrack);
- if (copy.PropagateToDCA(vertex, field, 100., b, bCov) )
- {
- vDCApropxy = b[0];
- vDCApropz = b[1];
- }
- else
- {
- vDCApropxy = -999999;
- vDCApropz = -999999;
- }
+ // kinematic cuts
+ if( vPt > fPtMax || vPt < fPtMin )
+ continue;
+ if( vEta > fEtaMax || vEta < fEtaMin )
+ continue;
+
+ // if not constrained track the position is stored (primary vertex to be subtracted)
+ aodTrack->GetXYZ(pos);
+ vDCAglobalx = pos[0] - v[0];
+ vDCAglobaly = pos[1] - v[1];
+ vDCAglobalz = pos[2] - v[2];
+ // fill for separately for positive and negative charges
+ Int_t iCharge = -1;
+ // positive
+ if(aodTrack->Charge() > 0)
+ iCharge = 0;
+ else if(aodTrack->Charge() < 0)
+ iCharge = 1;
+ else{
+ AliError("Charge==0?");
+ iCharge = -1;
+ }
+
+
// AOD track cuts
- for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
- fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
+ if(iCharge > -1){
+ for(Int_t iTrackBit = 0; iTrackBit < gBitMax-1; iTrackBit++){
+ fHistTrackStats->Fill(gCentrality,iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
+
+ if(aodTrack->TestFilterBit(1<<iTrackBit)){
+ fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
+ fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
+ fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
+ fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
+ }
+ }//bit loop
- if(aodTrack->TestFilterBit(1<<iTrackBit)){
- fHistKinematics[iTrackBit]->Fill(vEta,vPhi,vPt);
- fHistDCA[iTrackBit]->Fill(vDCAxy,vDCAz);
- fHistDCAprop[iTrackBit]->Fill(vDCApropxy,vDCApropz);
- fHistChiClus[iTrackBit]->Fill(vChi2,vClus);
- }
+ // fill all tracks in last bit
+ fHistTrackStats->Fill(gCentrality,gBitMax-1,1);
+ fHistKinematics[iCharge][gBitMax-1]->Fill(vEta,vPhi,vPt);
+ fHistDCAconstrained[iCharge][gBitMax-1]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
+ fHistDCAglobal[iCharge][gBitMax-1]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
+ fHistChiClus[iCharge][gBitMax-1]->Fill(vChi2,vClus);
- }//bit loop
+ }//charge positive or negative
}//track loop
+
+ return;
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskAODFilterBitQA::GetAcceptedHFVertexingTracks(AliVEvent *event, Double_t gCentrality){
+ // Checks track cuts (filter bits)
+ // from daughters of HF candidates
+ // Fills QA histograms
+
+ Double_t vEta;
+ Double_t vPhi;
+ Double_t vPt;
+ Double_t vDCAconstrainedxy;
+ Double_t vDCAconstrainedz;
+ Double_t vDCAglobalx;
+ Double_t vDCAglobaly;
+ Double_t vDCAglobalz;
+ Double_t vChi2;
+ Double_t vClus;
+
+ Double_t pos[3];
+ Double_t v[3];
+
+ Int_t IDList[1000];
+ Int_t IDListLength = 0;
+ for(Int_t i = 0; i < 1000; i++){
+ IDList[i] = -5;
+ }
+
+ const AliVVertex *vertex = event->GetPrimaryVertex();
+ vertex->GetXYZ(v);
+
+ // =================================================================================
+ // HF part (taken from AliAnalysisTaskSEDmesonsFilterCJ)
+
+ TClonesArray *arrayDStartoD0pi = (TClonesArray*)event->GetList()->FindObject(fHFBranchName.Data());
+
+ if (!arrayDStartoD0pi) {
+ AliInfo(Form("Could not find array %s, skipping the event",fHFBranchName.Data()));
+ return;
+ } else {
+ AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
+ }
+ // loop over tracks
+ const Int_t NVertices = arrayDStartoD0pi->GetEntriesFast();
+
+ for (Int_t iVertex = 0; iVertex < NVertices; ++iVertex) {
+
+ AliAODRecoDecay *HFvertex = static_cast<AliAODRecoDecay*>(arrayDStartoD0pi->At(iVertex));
+
+ // Loop over tracks (daughters of D candidates)
+ for (Int_t iProng = 0; iProng<HFvertex->GetNProngs(); iProng++) {
+
+ AliAODTrack *daughter = (AliAODTrack*)HFvertex->GetDaughter(iProng);
+ if (!daughter){
+ AliError(Form("Could not receive track %d %d", iVertex, iProng));
+ continue;
+ }
+
+
+ // get MC information (if available)
+ if(fArrayMC && fillOnlySecondaries){
+
+ Int_t label = daughter->GetLabel();
+ AliAODMCParticle *mcTrack = (AliAODMCParticle *)fArrayMC->At(TMath::Abs(label));
+
+ if(mcTrack->IsPhysicalPrimary())
+ continue;
+ }
+
+ // track parameters
+ vEta = daughter->Eta();
+ vPhi = daughter->Phi();// * TMath::RadToDeg();
+ vPt = daughter->Pt();
+ vDCAconstrainedxy = daughter->DCA();
+ vDCAconstrainedz = daughter->ZAtDCA();
+ vChi2 = daughter->Chi2perNDF();
+ vClus = daughter->GetTPCNcls();
+
+ // kinematic cuts
+ if( vPt > fPtMax || vPt < fPtMin )
+ continue;
+ if( vEta > fEtaMax || vEta < fEtaMin )
+ continue;
+
+ // avoid double counting (can be optimized)
+ Bool_t doubleCount = kFALSE;
+ Int_t daughterID = daughter->GetID();
+ for(Int_t idx = 0; idx < IDListLength; idx++){
+ if(IDList[idx]==daughterID){
+ doubleCount = kTRUE;
+ break;
+ }
+ }
+ if(!doubleCount){
+ IDList[IDListLength] = daughterID;
+ IDListLength++;
+ }
+ else{
+ continue;
+ }
+
+
+ // if not constrained track the position is stored (primary vertex to be subtracted)
+ daughter->GetXYZ(pos);
+ vDCAglobalx = pos[0] - v[0];
+ vDCAglobaly = pos[1] - v[1];
+ vDCAglobalz = pos[2] - v[2];
+
+ // fill for separately for positive and negative charges
+ Int_t iCharge = -1;
+ // positive
+ if(daughter->Charge() > 0)
+ iCharge = 0;
+ else if(daughter->Charge() < 0)
+ iCharge = 1;
+ else{
+ AliError("Charge==0?");
+ iCharge = -1;
+ }
+
+
+ // AOD track cuts
+ if(iCharge > -1){
+
+ // if some filter bits should be ignored, skip them here
+ if(fBitIgnore1 > -1 && daughter->TestFilterBit(1<<fBitIgnore1))
+ continue;
+ if(fBitIgnore2 > -1 && daughter->TestFilterBit(1<<fBitIgnore2))
+ continue;
+
+ for(Int_t iTrackBit = 0; iTrackBit < gBitMax; iTrackBit++){
+ fHistTrackStats->Fill(gCentrality,iTrackBit,daughter->TestFilterBit(1<<iTrackBit));
+
+ if(daughter->TestFilterBit(1<<iTrackBit)){
+ fHistKinematics[iCharge][iTrackBit]->Fill(vEta,vPhi,vPt);
+ fHistDCAconstrained[iCharge][iTrackBit]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
+ fHistDCAglobal[iCharge][iTrackBit]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
+ fHistChiClus[iCharge][iTrackBit]->Fill(vChi2,vClus);
+ }
+ }//bit loop
+
+ // fill all tracks in last bit
+ fHistTrackStats->Fill(gCentrality,gBitMax-1,1);
+ fHistKinematics[iCharge][gBitMax-1]->Fill(vEta,vPhi,vPt);
+ fHistDCAconstrained[iCharge][gBitMax-1]->Fill(vDCAconstrainedxy,vDCAconstrainedz);
+ fHistDCAglobal[iCharge][gBitMax-1]->Fill(vDCAglobalx,vDCAglobaly,vDCAglobalz);
+ fHistChiClus[iCharge][gBitMax-1]->Fill(vChi2,vClus);
+
+ }//charge positive or negative
+
+ }//prong loop
+ }//HF vertex loop
+
return;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff