/* $Id: AliAnalysisTaskESDfilter.cxx 24535 2008-03-16 22:43:30Z fca $ */
#include <TChain.h>
-#include <TFile.h>
+#include <TTree.h>
+#include <TList.h>
#include <TArrayI.h>
+#include <TRandom.h>
+#include <TParticle.h>
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
+#include "AliStack.h"
#include "AliAODEvent.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
#include "AliESDInputHandler.h"
#include "AliAODHandler.h"
+#include "AliAODMCParticle.h"
#include "AliAnalysisFilter.h"
-#include "AliESDtrack.h"
#include "AliESDMuonTrack.h"
#include "AliESDVertex.h"
#include "AliESDv0.h"
AliAnalysisTaskSE(),
fTrackFilter(0x0),
fKinkFilter(0x0),
- fV0Filter(0x0)
+ fV0Filter(0x0),
+ fCascadeFilter(0x0),
+ fHighPthreshold(0),
+ fPtshape(0x0)
{
// Default constructor
}
AliAnalysisTaskSE(name),
fTrackFilter(0x0),
fKinkFilter(0x0),
- fV0Filter(0x0)
+ fV0Filter(0x0),
+ fCascadeFilter(0x0),
+ fHighPthreshold(0),
+ fPtshape(0x0)
{
// Constructor
}
Long64_t ientry = Entry();
if (fDebug > 0) printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
+ if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
+ if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
ConvertESDtoAOD();
}
void AliAnalysisTaskESDfilter::ConvertESDtoAOD() {
// ESD Filter analysis task executed for each event
+
AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
- AliESD* old = esd->GetAliESDOld();
-
+ AliESD* old = esd->GetAliESDOld();
+
+ // Fetch Stack for debuggging if available
+ AliStack *pStack = 0;
+ AliMCEventHandler *mcH = 0;
+ if(MCEvent()){
+ pStack = MCEvent()->Stack();
+ mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
+ }
// set arrays and pointers
Float_t posF[3];
Double_t pos[3];
Double_t p[3];
- Double_t p_pos[3];
- Double_t p_neg[3];
+ Double_t momPos[3];
+ Double_t momNeg[3];
+ Double_t momBach[3];
+ Double_t momPosAtV0vtx[3];
+ Double_t momNegAtV0vtx[3];
+ Double_t momBachAtCascadeVtx[3];
Double_t covVtx[6];
Double_t covTr[21];
Double_t pid[10];
for (Int_t i = 0; i < 21; i++) covTr [i] = 0.;
- // loop over events and fill them
-
- // Multiplicity information needed by the header (to be revised!)
+ // loop over events and fill them
+
+ // Multiplicity information needed by the header (to be revised!)
Int_t nTracks = esd->GetNumberOfTracks();
+ // if (fDebug > 0) printf("-------------------Bo: Number of ESD tracks %d \n",nTracks);
+
Int_t nPosTracks = 0;
// for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
// if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
// Update the header
AliAODHeader* header = AODEvent()->GetHeader();
+
header->SetRunNumber(esd->GetRunNumber());
if (old) {
header->SetBunchCrossNumber(0);
header->SetZDCN2Energy(esd->GetZDCN2Energy());
header->SetZDCP2Energy(esd->GetZDCP2Energy());
header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
+ Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
+ Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
+ header->SetDiamond(diamxy,diamcov);
+//
//
-//
Int_t nV0s = esd->GetNumberOfV0s();
Int_t nCascades = esd->GetNumberOfCascades();
Int_t nKinks = esd->GetNumberOfKinks();
- Int_t nVertices = nV0s + 2*nCascades /*could lead to two vertices, one V0 and the Xi */+ nKinks + 1 /* = prim. vtx*/;
+ Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
Int_t nJets = 0;
Int_t nCaloClus = esd->GetNumberOfCaloClusters();
Int_t nFmdClus = 0;
if (fDebug > 0)
printf(" NV0=%d NCASCADES=%d NKINKS=%d\n", nV0s, nCascades, nKinks);
+
+ AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
- AODEvent()->ResetStd(nTracks, nVertices, nV0s+nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
- AliAODTrack *aodTrack = 0x0;
+ AliAODTrack *aodTrack = 0x0;
+ AliAODPid *detpid = 0x0;
+ Double_t timezero = 0; //TO BE FIXED
+ AliAODVertex *vV0FromCascade = 0x0;
+ AliAODv0 *aodV0 = 0x0;
+ AliAODcascade *aodCascade = 0x0;
+ // RefArray to store the mapping between esd track number and newly created AOD-Track
+ TRefArray *aodTrackRefs = NULL;
+ if (nTracks > 0) aodTrackRefs = new TRefArray(nTracks);
+
+ // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0
+ TRefArray *aodV0VtxRefs = NULL;
+ if (nV0s > 0) aodV0VtxRefs = new TRefArray(nV0s);
+
+ // RefArray to store the mapping between esd V0 number and newly created AOD-V0
+ TRefArray *aodV0Refs = NULL;
+ if (nV0s > 0) aodV0Refs = new TRefArray(nV0s);
+
+
+
+
+
// Array to take into account the tracks already added to the AOD
Bool_t * usedTrack = NULL;
if (nTracks>0) {
usedTrack = new Bool_t[nTracks];
for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
}
- // Array to take into account the V0s already added to the AOD
+ // Array to take into account the V0s already added to the AOD (V0 within cascades)
Bool_t * usedV0 = NULL;
if (nV0s>0) {
usedV0 = new Bool_t[nV0s];
TClonesArray &tracks = *(AODEvent()->GetTracks());
Int_t jTracks=0;
+ // Access to the AOD container of Cascades
+ TClonesArray &cascades = *(AODEvent()->GetCascades());
+ Int_t jCascades=0;
+
// Access to the AOD container of V0s
- TClonesArray &V0s = *(AODEvent()->GetV0s());
+ TClonesArray &v0s = *(AODEvent()->GetV0s());
Int_t jV0s=0;
// Add primary vertex. The primary tracks will be defined
vtx->GetCovMatrix(covVtx); //covariance matrix
AliAODVertex * primary = new(vertices[jVertices++])
- AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, AliAODVertex::kPrimary);
+ AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
+ primary->SetName(vtx->GetName());
+ primary->SetTitle(vtx->GetTitle());
+
if (fDebug > 0) primary->Print();
// Create vertices starting from the most complex objects
Double_t chi2 = 0.;
- // Cascades
+ // Cascades (Modified by A.Maire - February 2009)
for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
- AliESDcascade *cascade = esd->GetCascade(nCascade);
-
- cascade->GetXYZ(pos[0], pos[1], pos[2]);
- if (!old) {
- chi2 = cascade->GetChi2Xi(); // = chi2/NDF since NDF = 2*2-3
- cascade->GetPosCovXi(covVtx);
- } else {
- chi2 = -999.;
- for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
- }
- // Add the cascade vertex
- AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
- covVtx,
- chi2,
- primary,
- nCascade,
- AliAODVertex::kCascade);
-
- primary->AddDaughter(vcascade);
+ if (fDebug > 1)
+ printf("\n ******** Cascade number : %d/%d *********\n", nCascade, nCascades);
+
+ // 0- Preparation
+ //
+ AliESDcascade *esdCascade = esd->GetCascade(nCascade);
+ Int_t idxPosFromV0Dghter = esdCascade->GetPindex();
+ Int_t idxNegFromV0Dghter = esdCascade->GetNindex();
+ Int_t idxBachFromCascade = esdCascade->GetBindex();
- // Add the V0 from the cascade. The ESD class have to be optimized...
- // Now we have to search for the corresponding Vo in the list of V0s
- // using the indeces of the positive and negative tracks
+ AliESDtrack *esdCascadePos = esd->GetTrack( idxPosFromV0Dghter);
+ AliESDtrack *esdCascadeNeg = esd->GetTrack( idxNegFromV0Dghter);
+ AliESDtrack *esdCascadeBach = esd->GetTrack( idxBachFromCascade);
+
+ // Identification of the V0 within the esdCascade (via both daughter track indices)
+ AliESDv0 * currentV0 = 0x0;
+ Int_t idxV0FromCascade = -1;
- Int_t posFromV0 = cascade->GetPindex();
- Int_t negFromV0 = cascade->GetNindex();
+ for (Int_t iV0=0; iV0<nV0s; ++iV0) {
-
- AliESDv0 * v0 = 0x0;
- Int_t indV0 = -1;
+ currentV0 = esd->GetV0(iV0);
+ Int_t posCurrentV0 = currentV0->GetPindex();
+ Int_t negCurrentV0 = currentV0->GetNindex();
- for (Int_t iV0=0; iV0<nV0s; ++iV0) {
-
- v0 = esd->GetV0(iV0);
- Int_t posV0 = v0->GetPindex();
- Int_t negV0 = v0->GetNindex();
-
- if (posV0==posFromV0 && negV0==negFromV0) {
- indV0 = iV0;
+ if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {
+ idxV0FromCascade = iV0;
break;
- }
+ }
}
+
+ if(idxV0FromCascade < 0){
+ printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");
+ continue;
+ }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"
+
+ if (fDebug > 1)
+ printf("Cascade %d - V0fromCascade ind : %d/%d\n", nCascade, idxV0FromCascade, nV0s);
+
+ AliESDv0 *esdV0FromCascade = esd->GetV0(idxV0FromCascade);
- AliAODVertex * vV0FromCascade = 0x0;
+
+ // 1 - Cascade selection
- if (indV0>-1 && !usedV0[indV0] ) {
-
- // the V0 exists in the array of V0s and is not used
-
- usedV0[indV0] = kTRUE;
-
- v0->GetXYZ(pos[0], pos[1], pos[2]);
- if (!old) {
- chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
- v0->GetPosCov(covVtx);
- } else {
- chi2 = -999.;
- for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
- }
+ // AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd->GetPrimaryVertex()));
+ // TList cascadeObjects;
+ // cascadeObjects.AddAt(esdV0FromCascade, 0);
+ // cascadeObjects.AddAt(esdCascadePos, 1);
+ // cascadeObjects.AddAt(esdCascadeNeg, 2);
+ // cascadeObjects.AddAt(esdCascade, 3);
+ // cascadeObjects.AddAt(esdCascadeBach, 4);
+ // cascadeObjects.AddAt(esdPrimVtx, 5);
+ //
+ // UInt_t selectCascade = 0;
+ // if (fCascadeFilter) {
+ // // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects);
+ // // FIXME AliESDCascadeCuts to be implemented ...
+ //
+ // // Here we may encounter a moot point at the V0 level
+ // // between the cascade selections and the V0 ones :
+ // // the V0 selected along with the cascade (secondary V0) may
+ // // usually be removed from the dedicated V0 selections (prim V0) ...
+ // // -> To be discussed !
+ //
+ // // this is a little awkward but otherwise the
+ // // list wants to access the pointer (delete it)
+ // // again when going out of scope
+ // delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
+ // esdPrimVtx = 0;
+ // if (!selectCascade)
+ // continue;
+ // }
+ // else{
+ // delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
+ // esdPrimVtx = 0;
+ // }
- vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
- covVtx,
- chi2,
- vcascade,
- indV0,
- AliAODVertex::kV0);
- } else {
-
- // the V0 doesn't exist in the array of V0s or was used
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
-// << " The V0 " << indV0
-// << " doesn't exist in the array of V0s or was used!" << endl;
-
- cascade->GetXYZ(pos[0], pos[1], pos[2]);
-
- if (!old) {
- chi2 = v0->GetChi2V0();
- cascade->GetPosCov(covVtx);
- } else {
- chi2 = -999.;
- for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
- }
+ // 2 - Add the cascade vertex
+
+ esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);
+ esdCascade->GetPosCovXi(covVtx);
+ chi2 = esdCascade->GetChi2Xi();
+
+ AliAODVertex *vCascade = new(vertices[jVertices++]) AliAODVertex( pos,
+ covVtx,
+ chi2, // FIXME = Chi2/NDF will be needed
+ primary,
+ nCascade, // id
+ AliAODVertex::kCascade);
+ primary->AddDaughter(vCascade);
- vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
- covVtx,
- chi2, // = chi2/NDF since NDF = 2*2-3 (AM)
- vcascade,
- indV0,
- AliAODVertex::kV0);
- vcascade->AddDaughter(vV0FromCascade);
+ if (fDebug > 2) {
+ printf("---- Cascade / Cascade Vertex (AOD) : \n");
+ vCascade->Print();
}
+
+
+ // 3 - Add the bachelor track from the cascade
- // Add the positive tracks from the V0
-
- if (posFromV0>-1 && !usedTrack[posFromV0]) {
-
- usedTrack[posFromV0] = kTRUE;
-
- AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
- esdTrack->GetPxPyPz(p_pos);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
- vV0FromCascade->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p_pos,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- vV0FromCascade,
- kTRUE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kSecondary)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
+ if (!usedTrack[idxBachFromCascade]) {
+
+ esdCascadeBach->GetPxPyPz(momBach);
+ esdCascadeBach->GetXYZ(pos);
+ esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);
+ esdCascadeBach->GetESDpid(pid);
+
+ usedTrack[idxBachFromCascade] = kTRUE;
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);
+ if(mcH)mcH->SelectParticle(esdCascadeBach->GetLabel());
+ aodTrack = new(tracks[jTracks++]) AliAODTrack(esdCascadeBach->GetID(),
+ esdCascadeBach->GetLabel(),
+ momBach,
+ kTRUE,
+ pos,
+ kFALSE, // Why kFALSE for "isDCA" ? FIXME
+ covTr,
+ (Short_t)esdCascadeBach->GetSign(),
+ esdCascadeBach->GetITSClusterMap(),
+ pid,
+ vCascade,
+ kTRUE, // usedForVtxFit = kFALSE ? FIXME
+ vtx->UsesTrack(esdCascadeBach->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
+ aodTrackRefs->AddAt(aodTrack,idxBachFromCascade);
+
+ if (esdCascadeBach->GetSign() > 0) nPosTracks++;
aodTrack->ConvertAliPIDtoAODPID();
- aodTrack->SetFlags(esdTrack->GetStatus());
+ aodTrack->SetFlags(esdCascadeBach->GetStatus());
+ SetAODPID(esdCascadeBach,aodTrack,detpid,timezero,esd->GetMagneticField());
}
else {
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
-// << " track " << posFromV0 << " has already been used!" << endl;
+ aodTrack = dynamic_cast<AliAODTrack*>( aodTrackRefs->At(idxBachFromCascade) );
+ }
+
+ vCascade->AddDaughter(aodTrack);
+
+ if (fDebug > 4) {
+ printf("---- Cascade / bach dghter : \n");
+ aodTrack->Print();
}
- // Add the negative tracks from the V0
+
+ // 4 - Add the V0 from the cascade.
+ // = V0vtx + both pos and neg daughter tracks + the aodV0 itself
+ //
+
+ if ( !usedV0[idxV0FromCascade] ) {
+ // 4.A - if VO structure hasn't been created yet
+
+ // 4.A.1 - Create the V0 vertex of the cascade
- if (negFromV0>-1 && !usedTrack[negFromV0]) {
-
- usedTrack[negFromV0] = kTRUE;
-
- AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
- esdTrack->GetPxPyPz(p_neg);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
- vV0FromCascade->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p_neg,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- vV0FromCascade,
- kTRUE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kSecondary)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
- aodTrack->ConvertAliPIDtoAODPID();
- aodTrack->SetFlags(esdTrack->GetStatus());
- }
- else {
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
-// << " track " << negFromV0 << " has already been used!" << endl;
- }
+ esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);
+ esdV0FromCascade->GetPosCov(covVtx);
+ chi2 = esdV0FromCascade->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3 ?
+
+ vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
+ covVtx,
+ chi2,
+ vCascade,
+ idxV0FromCascade, //id of ESDv0
+ AliAODVertex::kV0);
+ // Note:
+ // one V0 can be used by several cascades.
+ // So, one AOD V0 vtx can have several parent vtx.
+ // This is not directly allowed by AliAODvertex.
+ // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash
+ // but to a problem of consistency within AODEvent.
+ // -> See below paragraph 4.B, for the proposed treatment of such a case.
+
+ // Add the vV0FromCascade to the aodVOVtxRefs
+ aodV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);
+
+
+ // 4.A.2 - Add the positive tracks from the V0
+
+ esdCascadePos->GetPxPyPz(momPos);
+ esdCascadePos->GetXYZ(pos);
+ esdCascadePos->GetCovarianceXYZPxPyPz(covTr);
+ esdCascadePos->GetESDpid(pid);
+
+
+ if (!usedTrack[idxPosFromV0Dghter]) {
+ usedTrack[idxPosFromV0Dghter] = kTRUE;
+
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);
+ if(mcH) mcH->SelectParticle(esdCascadePos->GetLabel());
+ aodTrack = new(tracks[jTracks++]) AliAODTrack( esdCascadePos->GetID(),
+ esdCascadePos->GetLabel(),
+ momPos,
+ kTRUE,
+ pos,
+ kFALSE, // Why kFALSE for "isDCA" ? FIXME
+ covTr,
+ (Short_t)esdCascadePos->GetSign(),
+ esdCascadePos->GetITSClusterMap(),
+ pid,
+ vV0FromCascade,
+ kTRUE, // usedForVtxFit = kFALSE ? FIXME
+ vtx->UsesTrack(esdCascadePos->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
+ aodTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);
+
+ if (esdCascadePos->GetSign() > 0) nPosTracks++;
+ aodTrack->ConvertAliPIDtoAODPID();
+ aodTrack->SetFlags(esdCascadePos->GetStatus());
+ SetAODPID(esdCascadePos,aodTrack,detpid,timezero,esd->GetMagneticField());
+ }
+ else {
+ aodTrack = dynamic_cast<AliAODTrack*>(aodTrackRefs->At(idxPosFromV0Dghter));
+ }
+ vV0FromCascade->AddDaughter(aodTrack);
+
+
+ // 4.A.3 - Add the negative tracks from the V0
+
+ esdCascadeNeg->GetPxPyPz(momNeg);
+ esdCascadeNeg->GetXYZ(pos);
+ esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);
+ esdCascadeNeg->GetESDpid(pid);
+
+
+ if (!usedTrack[idxNegFromV0Dghter]) {
+ usedTrack[idxNegFromV0Dghter] = kTRUE;
+
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);
+ if(mcH)mcH->SelectParticle(esdCascadeNeg->GetLabel());
+ aodTrack = new(tracks[jTracks++]) AliAODTrack( esdCascadeNeg->GetID(),
+ esdCascadeNeg->GetLabel(),
+ momNeg,
+ kTRUE,
+ pos,
+ kFALSE, // Why kFALSE for "isDCA" ? FIXME
+ covTr,
+ (Short_t)esdCascadeNeg->GetSign(),
+ esdCascadeNeg->GetITSClusterMap(),
+ pid,
+ vV0FromCascade,
+ kTRUE, // usedForVtxFit = kFALSE ? FIXME
+ vtx->UsesTrack(esdCascadeNeg->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
+
+ aodTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);
+
+ if (esdCascadeNeg->GetSign() > 0) nPosTracks++;
+ aodTrack->ConvertAliPIDtoAODPID();
+ aodTrack->SetFlags(esdCascadeNeg->GetStatus());
+ SetAODPID(esdCascadeNeg,aodTrack,detpid,timezero,esd->GetMagneticField());
+ }
+ else {
+ aodTrack = dynamic_cast<AliAODTrack*>(aodTrackRefs->At(idxNegFromV0Dghter));
+ }
- // add it to the V0 array as well
- Double_t d0[2] = { -999., -99.};
- // counting is probably wrong
- new(V0s[jV0s++]) AliAODv0(vV0FromCascade, -999., -99., p_pos, p_neg, d0); // to be refined
+ vV0FromCascade->AddDaughter(aodTrack);
+
+
+ // 4.A.4 - Add the V0 from cascade to the V0 array
- // Add the bachelor track from the cascade
+ Double_t dcaV0Daughters = esdV0FromCascade->GetDcaV0Daughters();
+ Double_t dcaV0ToPrimVertex = esdV0FromCascade->GetD( esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetPrimaryVertex()->GetZ() );
+ esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] );
+ esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] );
- Int_t bachelor = cascade->GetBindex();
+ Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
+ dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD( esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetMagneticField()) );
+ dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD( esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetMagneticField()) );
+
+ aodV0 = new(v0s[jV0s++]) AliAODv0( vV0FromCascade,
+ dcaV0Daughters,
+ dcaV0ToPrimVertex,
+ momPosAtV0vtx,
+ momNegAtV0vtx,
+ dcaDaughterToPrimVertex);
+ // set the aod v0 on-the-fly status
+ aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());
+
+ // Add the aodV0 to the aodVORefs
+ aodV0Refs->AddAt(aodV0,idxV0FromCascade);
+
+ usedV0[idxV0FromCascade] = kTRUE;
+
+ } else {
+ // 4.B - if V0 structure already used
+
+ // Note :
+ // one V0 can be used by several cascades (frequent in PbPb evts) :
+ // same V0 which used but attached to different bachelor tracks
+ // -> aodVORefs and aodV0VtxRefs are needed.
+ // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.
- if(bachelor>-1 && !usedTrack[bachelor]) {
-
- usedTrack[bachelor] = kTRUE;
-
- AliESDtrack *esdTrack = esd->GetTrack(bachelor);
- esdTrack->GetPxPyPz(p);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
- vcascade->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- vcascade,
- kTRUE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kSecondary)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
- aodTrack->ConvertAliPIDtoAODPID();
- aodTrack->SetFlags(esdTrack->GetStatus());
+ vV0FromCascade = dynamic_cast<AliAODVertex*>( aodV0VtxRefs->At(idxV0FromCascade) );
+ aodV0 = dynamic_cast<AliAODv0*> ( aodV0Refs ->At(idxV0FromCascade) );
+
+ // - Treatment of the parent for such a "re-used" V0 :
+ // Insert the cascade that reuses the V0 vertex in the lineage chain
+ // Before : vV0 -> vCascade1 -> vPrimary
+ // - Hyp : cascade2 uses the same V0 as cascade1
+ // After : vV0 -> vCascade2 -> vCascade1 -> vPrimary
+
+ AliAODVertex *vCascadePreviousParent = dynamic_cast<AliAODVertex*> (vV0FromCascade->GetParent());
+ vV0FromCascade->SetParent(vCascade);
+ vCascade ->SetParent(vCascadePreviousParent);
+
+ if(fDebug > 2)
+ printf("---- Cascade / Lineage insertion\n"
+ "Parent of V0 vtx = Cascade vtx %p\n"
+ "Parent of the cascade vtx = Cascade vtx %p\n"
+ "Parent of the parent cascade vtx = Cascade vtx %p\n",
+ static_cast<void*> (vV0FromCascade->GetParent()),
+ static_cast<void*> (vCascade->GetParent()),
+ static_cast<void*> (vCascadePreviousParent->GetParent()) );
+
+ }// end if V0 structure already used
+
+ if (fDebug > 2) {
+ printf("---- Cascade / V0 vertex: \n");
+ vV0FromCascade->Print();
}
- else {
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
-// << " track " << bachelor << " has already been used!" << endl;
+
+ if (fDebug > 4) {
+ printf("---- Cascade / pos dghter : \n");
+ aodTrack->Print();
+ printf("---- Cascade / neg dghter : \n");
+ aodTrack->Print();
+ printf("---- Cascade / aodV0 : \n");
+ aodV0->Print();
}
+
+ // In any case (used V0 or not), add the V0 vertex to the cascade one.
+ vCascade->AddDaughter(vV0FromCascade);
- // Add the primary track of the cascade (if any)
+
+ // 5 - Add the primary track of the cascade (if any)
+
+
+ // 6 - Add the cascade to the AOD array of cascades
+
+ Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetMagneticField()) );
+
+ esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);
+
+ aodCascade = new(cascades[jCascades++]) AliAODcascade( vCascade,
+ esdCascade->Charge(),
+ esdCascade->GetDcaXiDaughters(),
+ -999.,
+ // DCAXiToPrimVtx -> needs to be calculated ----|
+ // doesn't exist at ESD level;
+ // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)
+ dcaBachToPrimVertexXY,
+ momBachAtCascadeVtx,
+ *aodV0);
+ if (fDebug > 3) {
+ printf("---- Cascade / AOD cascade : \n\n");
+ aodCascade->PrintXi(primary->GetX(), primary->GetY(), primary->GetZ());
+ }
+
} // end of the loop on cascades
-
+
+ cascades.Expand(jCascades);
+
+
+ //
// V0s
+ //
for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
- if (usedV0[nV0]) continue; // skip if aready added to the AOD
+ if (usedV0[nV0]) continue; // skip if already added to the AOD
AliESDv0 *v0 = esd->GetV0(nV0);
+ Int_t posFromV0 = v0->GetPindex();
+ Int_t negFromV0 = v0->GetNindex();
+ // V0 selection
+ //
+ AliESDVertex *esdVtx = new AliESDVertex(*(esd->GetPrimaryVertex()));
+ AliESDtrack *esdV0Pos = esd->GetTrack(posFromV0);
+ AliESDtrack *esdV0Neg = esd->GetTrack(negFromV0);
+ TList v0objects;
+ v0objects.AddAt(v0, 0);
+ v0objects.AddAt(esdV0Pos, 1);
+ v0objects.AddAt(esdV0Neg, 2);
+ v0objects.AddAt(esdVtx, 3);
+ UInt_t selectV0 = 0;
+ if (fV0Filter) {
+ selectV0 = fV0Filter->IsSelected(&v0objects);
+ // this is a little awkward but otherwise the
+ // list wants to access the pointer (delete it)
+ // again when going out of scope
+ delete v0objects.RemoveAt(3); // esdVtx created via copy construct
+ esdVtx = 0;
+ if (!selectV0)
+ continue;
+ }
+ else{
+ delete v0objects.RemoveAt(3); // esdVtx created via copy construct
+ esdVtx = 0;
+ }
+
v0->GetXYZ(pos[0], pos[1], pos[2]);
if (!old) {
AliAODVertex::kV0);
primary->AddDaughter(vV0);
- Int_t posFromV0 = v0->GetPindex();
- Int_t negFromV0 = v0->GetNindex();
-
+
// Add the positive tracks from the V0
- if (posFromV0>-1 && !usedTrack[posFromV0]) {
-
+ esdV0Pos->GetPxPyPz(momPos);
+ esdV0Pos->GetXYZ(pos);
+ esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
+ esdV0Pos->GetESDpid(pid);
+
+ if (!usedTrack[posFromV0]) {
usedTrack[posFromV0] = kTRUE;
-
- AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
- esdTrack->GetPxPyPz(p_pos);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
- vV0->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p_pos,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- vV0,
- kTRUE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kSecondary)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);
+ if(mcH)mcH->SelectParticle(esdV0Pos->GetLabel());
+ aodTrack = new(tracks[jTracks++]) AliAODTrack(esdV0Pos->GetID(),
+ esdV0Pos->GetLabel(),
+ momPos,
+ kTRUE,
+ pos,
+ kFALSE,
+ covTr,
+ (Short_t)esdV0Pos->GetSign(),
+ esdV0Pos->GetITSClusterMap(),
+ pid,
+ vV0,
+ kTRUE, // check if this is right
+ vtx->UsesTrack(esdV0Pos->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
+ aodTrackRefs->AddAt(aodTrack,posFromV0);
+ // if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
+ if (esdV0Pos->GetSign() > 0) nPosTracks++;
aodTrack->ConvertAliPIDtoAODPID();
- aodTrack->SetFlags(esdTrack->GetStatus());
+ aodTrack->SetFlags(esdV0Pos->GetStatus());
+ SetAODPID(esdV0Pos,aodTrack,detpid,timezero,esd->GetMagneticField());
}
else {
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " V0 " << nV0
-// << " track " << posFromV0 << " has already been used!" << endl;
+ aodTrack = dynamic_cast<AliAODTrack*>(aodTrackRefs->At(posFromV0));
+ // if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
}
-
+ vV0->AddDaughter(aodTrack);
+
// Add the negative tracks from the V0
- if (negFromV0>-1 && !usedTrack[negFromV0]) {
-
+ esdV0Neg->GetPxPyPz(momNeg);
+ esdV0Neg->GetXYZ(pos);
+ esdV0Neg->GetCovarianceXYZPxPyPz(covTr);
+ esdV0Neg->GetESDpid(pid);
+
+ if (!usedTrack[negFromV0]) {
usedTrack[negFromV0] = kTRUE;
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);
+ if(mcH)mcH->SelectParticle(esdV0Neg->GetLabel());
+ aodTrack = new(tracks[jTracks++]) AliAODTrack(esdV0Neg->GetID(),
+ esdV0Neg->GetLabel(),
+ momNeg,
+ kTRUE,
+ pos,
+ kFALSE,
+ covTr,
+ (Short_t)esdV0Neg->GetSign(),
+ esdV0Neg->GetITSClusterMap(),
+ pid,
+ vV0,
+ kTRUE, // check if this is right
+ vtx->UsesTrack(esdV0Neg->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
- AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
- esdTrack->GetPxPyPz(p_neg);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
- vV0->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p_neg,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- vV0,
- kTRUE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kSecondary)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
+ aodTrackRefs->AddAt(aodTrack,negFromV0);
+ // if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
+ if (esdV0Neg->GetSign() > 0) nPosTracks++;
aodTrack->ConvertAliPIDtoAODPID();
- aodTrack->SetFlags(esdTrack->GetStatus());
+ aodTrack->SetFlags(esdV0Neg->GetStatus());
+ SetAODPID(esdV0Neg,aodTrack,detpid,timezero,esd->GetMagneticField());
}
else {
-// cerr << "Error: event " << esd->GetEventNumberInFile() << " V0 " << nV0
-// << " track " << negFromV0 << " has already been used!" << endl;
+ aodTrack = dynamic_cast<AliAODTrack*>(aodTrackRefs->At(negFromV0));
+ // if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
}
+ vV0->AddDaughter(aodTrack);
+
+
+ // Add the V0 the V0 array as well
+
+ Double_t dcaV0Daughters = v0->GetDcaV0Daughters();
+ Double_t dcaV0ToPrimVertex = v0->GetD(esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetPrimaryVertex()->GetZ());
+ v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]);
+ v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]);
+
+ Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
+ dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD( esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetMagneticField()) );
+ dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD( esd->GetPrimaryVertex()->GetX(),
+ esd->GetPrimaryVertex()->GetY(),
+ esd->GetMagneticField()) );
+
+ aodV0 = new(v0s[jV0s++]) AliAODv0(vV0,
+ dcaV0Daughters,
+ dcaV0ToPrimVertex,
+ momPosAtV0vtx,
+ momNegAtV0vtx,
+ dcaDaughterToPrimVertex);
+
+ // set the aod v0 on-the-fly status
+ aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
+ }//End of loop on V0s
+
+ v0s.Expand(jV0s);
+
+ if (fDebug > 0) printf(" NAODCascades=%d / NAODV0s=%d\n", jCascades, jV0s);
+ // end of V0 parts
+
- // add it to the V0 array as well
- Double_t d0[2] = { 999., 99.};
- new(V0s[jV0s++]) AliAODv0(vV0, 999., 99., p_pos, p_neg, d0); // to be refined
- }
- V0s.Expand(jV0s);
- // end of the loop on V0s
-
// Kinks: it is a big mess the access to the information in the kinks
// The loop is on the tracks in order to find the mother and daugther of each kink
continue;
}
- // Add the mother track
+ // Add the mother track if it passed primary track selection cuts
AliAODTrack * mother = NULL;
+ UInt_t selectInfo = 0;
+ if (fTrackFilter) {
+ selectInfo = fTrackFilter->IsSelected(esd->GetTrack(imother));
+ if (!selectInfo) continue;
+ }
+
if (!usedTrack[imother]) {
usedTrack[imother] = kTRUE;
- AliESDtrack *esdTrack = esd->GetTrack(imother);
- esdTrack->GetPxPyPz(p);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
+ AliESDtrack *esdTrackM = esd->GetTrack(imother);
+ esdTrackM->GetPxPyPz(p);
+ esdTrackM->GetXYZ(pos);
+ esdTrackM->GetCovarianceXYZPxPyPz(covTr);
+ esdTrackM->GetESDpid(pid);
+ if(mcH)mcH->SelectParticle(esdTrackM->GetLabel());
mother =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
+ new(tracks[jTracks++]) AliAODTrack(esdTrackM->GetID(),
+ esdTrackM->GetLabel(),
p,
kTRUE,
pos,
kFALSE,
covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
+ (Short_t)esdTrackM->GetSign(),
+ esdTrackM->GetITSClusterMap(),
pid,
primary,
kTRUE, // check if this is right
- kTRUE, // check if this is right
- AliAODTrack::kPrimary);
- if (esdTrack->GetSign() > 0) nPosTracks++;
- mother->SetFlags(esdTrack->GetStatus());
+ vtx->UsesTrack(esdTrack->GetID()),
+ AliAODTrack::kPrimary,
+ selectInfo);
+ aodTrackRefs->AddAt(mother, imother);
+
+ if (esdTrackM->GetSign() > 0) nPosTracks++;
+ mother->SetFlags(esdTrackM->GetStatus());
mother->ConvertAliPIDtoAODPID();
primary->AddDaughter(mother);
mother->ConvertAliPIDtoAODPID();
+ SetAODPID(esdTrackM,mother,detpid,timezero,esd->GetMagneticField());
}
else {
// cerr << "Error: event " << esd->GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
usedTrack[idaughter] = kTRUE;
- AliESDtrack *esdTrack = esd->GetTrack(idaughter);
- esdTrack->GetPxPyPz(p);
- esdTrack->GetXYZ(pos);
- esdTrack->GetCovarianceXYZPxPyPz(covTr);
- esdTrack->GetESDpid(pid);
-
+ AliESDtrack *esdTrackD = esd->GetTrack(idaughter);
+ esdTrackD->GetPxPyPz(p);
+ esdTrackD->GetXYZ(pos);
+ esdTrackD->GetCovarianceXYZPxPyPz(covTr);
+ esdTrackD->GetESDpid(pid);
+ selectInfo = 0;
+ if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);
+ if(mcH)mcH->SelectParticle(esdTrackD->GetLabel());
daughter =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
+ new(tracks[jTracks++]) AliAODTrack(esdTrackD->GetID(),
+ esdTrackD->GetLabel(),
p,
kTRUE,
pos,
kFALSE,
covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
+ (Short_t)esdTrackD->GetSign(),
+ esdTrackD->GetITSClusterMap(),
pid,
vkink,
kTRUE, // check if this is right
- kTRUE, // check if this is right
- AliAODTrack::kPrimary);
- if (esdTrack->GetSign() > 0) nPosTracks++;
- daughter->SetFlags(esdTrack->GetStatus());
+ vtx->UsesTrack(esdTrack->GetID()),
+ AliAODTrack::kSecondary,
+ selectInfo);
+
+ aodTrackRefs->AddAt(daughter, idaughter);
+
+ if (esdTrackD->GetSign() > 0) nPosTracks++;
+ daughter->SetFlags(esdTrackD->GetStatus());
daughter->ConvertAliPIDtoAODPID();
vkink->AddDaughter(daughter);
daughter->ConvertAliPIDtoAODPID();
+ SetAODPID(esdTrackD,daughter,detpid,timezero,esd->GetMagneticField());
}
else {
// cerr << "Error: event " << esd->GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
// Track selection
if (fTrackFilter) {
selectInfo = fTrackFilter->IsSelected(esdTrack);
- if (!selectInfo) continue;
+ if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;
}
//
esdTrack->GetXYZ(pos);
esdTrack->GetCovarianceXYZPxPyPz(covTr);
esdTrack->GetESDpid(pid);
+ if(mcH)mcH->SelectParticle(esdTrack->GetLabel());
+ primary->AddDaughter(aodTrack =
+ new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
+ esdTrack->GetLabel(),
+ p,
+ kTRUE,
+ pos,
+ kFALSE,
+ covTr,
+ (Short_t)esdTrack->GetSign(),
+ esdTrack->GetITSClusterMap(),
+ pid,
+ primary,
+ kTRUE, // check if this is right
+ vtx->UsesTrack(esdTrack->GetID()),
+ AliAODTrack::kPrimary,
+ selectInfo)
+ );
+ aodTrackRefs->AddAt(aodTrack, nTrack);
- Float_t impactXY, impactZ;
-
- esdTrack->GetImpactParameters(impactXY,impactZ);
-
- if (impactXY<3) {
- // track inside the beam pipe
-
- primary->AddDaughter(aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- primary,
- kTRUE, // check if this is right
- kTRUE, // check if this is right
- AliAODTrack::kPrimary,
- selectInfo)
- );
- if (esdTrack->GetSign() > 0) nPosTracks++;
- aodTrack->SetFlags(esdTrack->GetStatus());
- aodTrack->ConvertAliPIDtoAODPID();
- }
- else {
- // outside the beam pipe: orphan track
- // Don't write them anymore!
- continue;
- /*
- aodTrack =
- new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
- esdTrack->GetLabel(),
- p,
- kTRUE,
- pos,
- kFALSE,
- covTr,
- (Short_t)esdTrack->GetSign(),
- esdTrack->GetITSClusterMap(),
- pid,
- NULL,
- kFALSE, // check if this is right
- kFALSE, // check if this is right
- AliAODTrack::kOrphan,
- selectInfo);
- if (esdTrack->GetSign() > 0) nPosTracks++;
- aodTrack->SetFlags(esdTrack->GetStatus());
- aodTrack->ConvertAliPIDtoAODPID();
- */
- }
+ if (esdTrack->GetSign() > 0) nPosTracks++;
+ aodTrack->SetFlags(esdTrack->GetStatus());
+ aodTrack->ConvertAliPIDtoAODPID();
+ SetAODPID(esdTrack,aodTrack,detpid,timezero,esd->GetMagneticField());
} // end of loop on tracks
// Update number of AOD tracks in header at the end of track loop (M.G.)
if (fDebug > 0)
printf(" NAODTRACKS=%d NPOS=%d NNEG=%d\n", jTracks, nPosTracks, jTracks - nPosTracks);
// Do not shrink the array of tracks - other filters may add to it (M.G)
-// tracks.Expand(jTracks); // remove 'empty slots' due to unwritten tracks
+ // tracks.Expand(jTracks); // remove 'empty slots' due to unwritten tracks
// Access to the AOD container of PMD clusters
TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
AliESDPmdTrack *pmdTrack = esd->GetPmdTrack(iPmd);
Int_t nLabel = 0;
Int_t *label = 0x0;
- Double_t pos[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ() };
- Double_t pid[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; // to be revised!
+ Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
+ Double_t pidPmd[13] = { 0.}; // to be revised!
// type not set!
// assoc cluster not set
- new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), pos, pid);
+ new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
}
// Access to the AOD container of clusters
Int_t nLabel = cluster->GetNLabels();
TArrayI* labels = cluster->GetLabels();
Int_t *label = 0;
- if (labels) label = (cluster->GetLabels())->GetArray();
+ if (labels){
+ label = (cluster->GetLabels())->GetArray();
+ for(int i = 0;i < labels->GetSize();++i){
+ if(mcH)mcH->SelectParticle(label[i]);
+ }
+ }
Float_t energy = cluster->E();
cluster->GetPosition(posF);
nLabel,
label,
energy,
- pos,
+ posF,
NULL,
ttype);
cluster->GetEmcCpvDistance(),
cluster->GetNExMax(),cluster->GetTOF()) ;
+ caloCluster->SetPIDFromESD(cluster->GetPid());
caloCluster->SetNCells(cluster->GetNCells());
caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
TArrayI* matchedT = cluster->GetTracksMatched();
- if (matchedT) {
+ if (nTracks>0 && matchedT && cluster->GetTrackMatched() >= 0) {
for (Int_t im = 0; im < matchedT->GetSize(); im++) {
- caloCluster->AddTrackMatched((esd->GetTrack(im)));
+ Int_t iESDtrack = matchedT->At(im);;
+ if (aodTrackRefs->At(iESDtrack) != 0) {
+ caloCluster->AddTrackMatched((AliAODTrack*)aodTrackRefs->At(iESDtrack));
+ }
}
}
-
+
}
caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
// end of loop on calo clusters
SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
- SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0), mult->GetLabel(n, 1));
+ if(mcH){
+ mcH->SelectParticle(mult->GetLabel(n, 0));
+ mcH->SelectParticle(mult->GetLabel(n, 1));
+ }
+ SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
}
}
} else {
//Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
}
- delete [] usedTrack;
- delete [] usedV0;
delete [] usedKink;
+ delete [] usedV0;
+ delete [] usedTrack;
+ delete aodV0Refs;
+ delete aodV0VtxRefs;
+ delete aodTrackRefs;
return;
}
+
+void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid, Double_t timezero, Double_t bfield)
+{
+ //
+ // Setter for the raw PID detector signals
+ //
+
+ if(esdtrack->Pt()>fHighPthreshold) {
+ detpid = new AliAODPid();
+ SetDetectorRawSignals(detpid,esdtrack,timezero, bfield);
+ aodtrack->SetDetPID(detpid);
+ } else {
+ if(fPtshape){
+ if(esdtrack->Pt()> fPtshape->GetXmin()){
+ Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
+ if(gRandom->Rndm(0)<1./y){
+ detpid = new AliAODPid();
+ SetDetectorRawSignals(detpid,esdtrack,timezero, bfield);
+ aodtrack->SetDetPID(detpid);
+ }//end rndm
+ }//end if p < pmin
+ }//end if p function
+ }// end else
+}
+
+void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track, Double_t timezero, Double_t bfield)
+{
+//
+//assignment of the detector signals (AliXXXesdPID inspired)
+//
+ if(!track){
+ AliInfo("no ESD track found. .....exiting");
+ return;
+ }
+
+ aodpid->SetITSsignal(track->GetITSsignal());
+ aodpid->SetTPCsignal(track->GetTPCsignal());
+ //n TRD planes = 6
+
+ Int_t nslices = track->GetNumberOfTRDslices()*6;
+ Double_t *trdslices = new Double_t[nslices];
+ for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
+ for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
+ }
+
+
+ aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);
+ Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
+ aodpid->SetIntegratedTimes(times);
+
+ aodpid->SetTOFsignal(track->GetTOFsignal()-timezero); // to be fixed
+ aodpid->SetHMPIDsignal(track->GetHMPIDsignal());
+
+ //Extrapolate track to EMCAL surface for AOD-level track-cluster matching
+ Double_t emcpos[3] = {0.,0.,0.};
+ Double_t emcmom[3] = {0.,0.,0.};
+ aodpid->SetEMCALPosition(emcpos);
+ aodpid->SetEMCALMomentum(emcmom);
+
+ AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();
+ if(!outerparam) return;
+
+ //To be replaced by call to AliEMCALGeoUtils when the class becomes available
+ Double_t radius = 441.0; //[cm] EMCAL radius +13cm
+
+ Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos);
+ Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom);
+ if(!(okpos && okmom)) return;
+
+ aodpid->SetEMCALPosition(emcpos);
+ aodpid->SetEMCALMomentum(emcmom);
+
+}
+
void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
{
// Terminate analysis
if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
}
+void AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){
+ if(!pStack)return;
+ label = TMath::Abs(label);
+ TParticle *part = pStack->Particle(label);
+ Printf("########################");
+ Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
+ part->Print();
+ TParticle* mother = part;
+ Int_t imo = part->GetFirstMother();
+ Int_t nprim = pStack->GetNprimary();
+ // while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
+ while((imo >= nprim)) {
+ mother = pStack->Particle(imo);
+ Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
+ mother->Print();
+ imo = mother->GetFirstMother();
+ }
+ Printf("########################");
+}