#include <TString.h>
#include <TMath.h>
+#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "AliAODVertex.h"
#include "AliVParticle.h"
#include "AliVertexerTracks.h"
#include "AliVVertex.h"
+#include "AliExternalTrackParam.h"
#include "AliHFEextraCuts.h"
ClassImp(AliHFEextraCuts)
-const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
+const Int_t AliHFEextraCuts::fgkNQAhistos = 10;
//______________________________________________________
AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
fMinNClustersTPCPID(0),
fClusterRatioTPC(0.),
fMinTrackletsTRD(0),
+ fMaxChi2TRD(5.0),
fMinNbITScls(0),
fTRDtrackletsExact(0),
fPixelITS(0),
+ fDriftITS(0),
fTPCclusterDef(0),
fTPCclusterRatioDef(0),
fFractionTPCShared(-1.0),
- fAbsHFEImpactParamNsigmaR(kTRUE),
+ fOppSideIPcut(kFALSE),
+ fTOFsignalDx(1.0),
+ fTOFsignalDz(1.0),
+ fMagField(-10000),
+ fAODFilterBit(0),
fCheck(kFALSE),
fQAlist(0x0) ,
fDebugLevel(0)
fMinNClustersTPCPID(c.fMinNClustersTPCPID),
fClusterRatioTPC(c.fClusterRatioTPC),
fMinTrackletsTRD(c.fMinTrackletsTRD),
+ fMaxChi2TRD(c.fMaxChi2TRD),
fMinNbITScls(c.fMinNbITScls),
fTRDtrackletsExact(c.fTRDtrackletsExact),
fPixelITS(c.fPixelITS),
+ fDriftITS(c.fDriftITS),
fTPCclusterDef(c.fTPCclusterDef),
fTPCclusterRatioDef(c.fTPCclusterRatioDef),
fFractionTPCShared(c.fFractionTPCShared),
- fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
+ fOppSideIPcut(c.fOppSideIPcut),
+ fTOFsignalDx(c.fTOFsignalDx),
+ fTOFsignalDz(c.fTOFsignalDz),
+ fMagField(c.fMagField),
+ fAODFilterBit(c.fAODFilterBit),
fCheck(c.fCheck),
fQAlist(0x0),
fDebugLevel(0)
fMinNClustersTPC = c.fMinNClustersTPC;
fMinNClustersTPCPID = c.fMinNClustersTPCPID;
fMinTrackletsTRD = c.fMinTrackletsTRD;
+ fMaxChi2TRD = c.fMaxChi2TRD;
fMinNbITScls = c.fMinNbITScls;
fTRDtrackletsExact = c.fTRDtrackletsExact;
fPixelITS = c.fPixelITS;
+ fDriftITS = c.fDriftITS;
fTPCclusterDef = c.fTPCclusterDef;
fTPCclusterRatioDef = c.fTPCclusterRatioDef;
fFractionTPCShared = c.fFractionTPCShared;
- fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
+ fOppSideIPcut = c.fOppSideIPcut;
+ fTOFsignalDx = c.fTOFsignalDx;
+ fTOFsignalDz = c.fTOFsignalDz;
+ fMagField = c.fMagField;
+ fAODFilterBit = c.fAODFilterBit;
fCheck = c.fCheck;
fDebugLevel = c.fDebugLevel;
memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
// QA histograms are filled before track selection and for
// selected tracks after track selection
//
- AliDebug(1, "Called");
+ AliDebug(1, Form("%s: Called", GetName()));
ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
// Apply cuts
- Float_t impactR, impactZ;
+ Float_t impactR = -999.;
+ Float_t impactZ = -999.;
Double_t hfeimpactR, hfeimpactnsigmaR;
Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
Double_t maximpactRcut;
GetImpactParameters(track, impactR, impactZ);
- if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
+ if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
// Protection for PbPb
GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
Double_t ratioTPC = GetTPCclusterRatio(track);
UChar_t trdTracklets;
trdTracklets = track->GetTRDntrackletsPID();
+ Float_t trdchi2=-999.;
+ trdchi2=GetTRDchi(track);
UChar_t itsPixel = track->GetITSClusterMap();
Int_t status1 = GetITSstatus(track, 0);
Int_t status2 = GetITSstatus(track, 1);
Bool_t statusL0 = CheckITSstatus(status1);
Bool_t statusL1 = CheckITSstatus(status2);
+ Double_t tofsignalDx = 0.0;
+ Double_t tofsignalDz = 0.0;
+ GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
+
if(TESTBIT(fRequirements, kTPCfractionShared)) {
// cut on max fraction of shared TPC clusters
if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
// cut on min. HFE Impact Parameter in Radial direction
if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
}
+ if(TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
+ // cut on min. HFE Impact Parameter in Radial direction, multiplied by particle charge
+ Double_t charge = (Double_t)track->Charge();
+
+ if(fMagField < 0)charge *= -1.;//the IP distribution side to be chosen depends on magnetic field polarization
+ if(fOppSideIPcut)charge*=-1.;//in case we choose to select electrons from the side of the photon peak
+ Double_t hfeimpactRcharge = hfeimpactR*charge;//switch selected side of the peak
+ if(hfeimpactRcharge >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamRcharge);
+ }
if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
// cut on max. HFE Impact Parameter n sigma in Radial direction
- if(fAbsHFEImpactParamNsigmaR) {
-// if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
- if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
+ // if(fAbsHFEImpactParamNsigmaR) {
+ // //if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
+ // if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
+ // SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+ // // printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+ // }
+ // }
+ // else {
+ if(hfeimpactnsigmaRcut>0){
+ if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
-// printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+ //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
}
}
- else {
- if(hfeimpactnsigmaRcut>0){
- if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
- SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
- //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
- }
- }
- else{
- if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
- SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
- //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
- }
+ else{
+ if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
+ SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+ //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
}
}
+ //}
}
if(TESTBIT(fRequirements, kClusterRatioTPC)){
// cut on min ratio of found TPC clusters vs findable TPC clusters
}
if(TESTBIT(fRequirements, kMinTrackletsTRD)){
// cut on minimum number of TRD tracklets
- AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
+ AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
if(fTRDtrackletsExact){
if(trdTracklets == fMinTrackletsTRD) {
SETBIT(survivedCut, kMinTrackletsTRD);
- AliDebug(1, "Track Selected");
+ AliDebug(1, Form("%s: Track Selected", GetName()));
}
}else{
if(trdTracklets >= fMinTrackletsTRD){
SETBIT(survivedCut, kMinTrackletsTRD);
- AliDebug(1, "Track Selected");
+ AliDebug(1, Form("%s: Track Selected", GetName()));
}
//printf("Min number of tracklets %d\n",fMinTrackletsTRD);
}
}
-
+
+ if(TESTBIT(fRequirements, kMaxTRDChi2)){
+ // cut on TRD chi2
+ AliDebug(1, Form("%s: Cut on TRD chi2: [%f|%f]\n", GetName(),fMaxChi2TRD, trdchi2));
+ if(trdchi2 < fMaxChi2TRD) {
+ SETBIT(survivedCut, kMaxTRDChi2);
+ AliDebug(1,Form("%s: survived %f\n",GetName(),trdchi2));
+ }
+ }
+
if(TESTBIT(fRequirements, kMinNbITScls)){
// cut on minimum number of ITS clusters
//printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
- AliDebug(1, Form("Min ITS clusters: [%d|%d]\n", fMinNbITScls, nclsITS));
+ AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
}
if(TESTBIT(fRequirements, kMinNClustersTPC)){
// cut on minimum number of TPC tracklets
//printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
- AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
+ AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
}
if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
- AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
+ AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
}
+ if(TESTBIT(fRequirements, kDriftITS)){
+ //printf("Require drift\n");
+ switch(fDriftITS){
+ case kFirstD:
+ if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
+ break;
+ default:
+ SETBIT(survivedCut, kDriftITS);
+ break;
+ }
+ }
if(TESTBIT(fRequirements, kPixelITS)){
// cut on ITS pixel layers
- AliDebug(1, "ITS cluster Map: ");
+ AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
//PrintBitMap(itsPixel);
switch(fPixelITS){
case kFirst:
SETBIT(survivedCut, kPixelITS);
break;
}
- AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
+ AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
}
if(TESTBIT(fRequirements, kTOFPID)){
if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
//else printf("Found kink mother\n");
}
+ if(TESTBIT(fRequirements, kTOFsignalDxy)){
+ // cut on TOF matching cluster
+ if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
+ }
+ if(TESTBIT(fRequirements, kITSpattern)){
+ // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
+ if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern);
+ }
+ if(TESTBIT(fRequirements, kAODFilterBit)){
+ AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+ if(aodtrack){
+ Int_t aodfilter = 1 << fAODFilterBit;
+ if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
+ }
+ }
if(fRequirements == survivedCut){
- //
+ //
// Track selected
//
- AliDebug(2, "Track Survived cuts\n");
+ AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
return kTRUE;
}
- AliDebug(2, "Track cut");
+ AliDebug(2, Form("%s: Track cut", GetName()));
if(IsQAOn()) FillCutCorrelation(survivedCut);
return kFALSE;
}
if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
}
if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
+
+ if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
+
+ if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
}
// //______________________________________________________
fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
histo1D->GetXaxis()->SetTitle("Impact Parameter");
histo1D->GetYaxis()->SetTitle("Number of Tracks");
- qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
+ qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
histo1D->GetYaxis()->SetTitle("Number of Tracks");
fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
histo1D->GetYaxis()->SetTitle("counts");
+ qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
+ fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
+ histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
+ histo1D->GetYaxis()->SetTitle("Number of Tracks");
+ qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
+ fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
+ histo1D->GetXaxis()->SetTitle("ITScls");
+ histo1D->GetYaxis()->SetTitle("Number of ITS cls");
+
}
// Add cut correlation
qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
}
//______________________________________________________
-Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
+Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
//
// Check ITS layer status
//
if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
Int_t det;
Float_t xloc, zloc;
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+ const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
}
return status;
//______________________________________________________
UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
- //
- // Get Number of findable clusters in the TPC
- //
- Int_t nClusters = 0; // in case no Information available consider all clusters findable
+ //
+ // Get Number of findable clusters in the TPC
+ //
+ Int_t nClusters = 0; // in case no Information available consider all clusters findable
TClass *type = track->IsA();
if(TESTBIT(fTPCclusterDef, kFoundIter1)){
- if(type == AliESDtrack::Class()){
+ if(type == AliESDtrack::Class()){
AliDebug(2, ("Using def kFoundIter1"));
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
nClusters = esdtrack->GetTPCNclsIter1();
} else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
AliDebug(2, ("Using def kCrossedRows"));
nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
- } else{
+ } else if(TESTBIT(fTPCclusterDef, kFound)){
AliDebug(2, ("Using def kFound"));
- nClusters = track->GetTPCNcls();
+ nClusters = track->GetTPCNcls();
+ }
+ else {
+ AliDebug(2, ("Using def kFoundAll"));
+ if(type == AliESDtrack::Class()){
+ AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+ const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
+ nClusters = clusterTPC.CountBits();
+ }
+ else {
+ AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+ const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
+ nClusters = clusterTPC.CountBits();
+ }
}
- return nClusters;
+ return nClusters;
}
//______________________________________________________
if(type == AliESDtrack::Class()){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
nClustersTPCShared = esdtrack->GetTPCnclsS();
- } else if(type == AliESDtrack::Class()){
+ } else if(type == AliAODTrack::Class()){
AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
const TBits &shared = aodtrack->GetTPCSharedMap();
nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
}
if (nClustersTPC!=0) {
- fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
+ fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
}
return fracClustersTPCShared;
}
AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
clusterRatio = 1.;
}
- } else {
+ } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
AliDebug(2, "Using ratio def kFoundOverFindable");
clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
}
+ else {
+ AliDebug(2, "Using ratio def kFoundAllOverFindable");
+ Int_t allclusters = 0;
+ if(type == AliESDtrack::Class()){
+ AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+ const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
+ allclusters = clusterTPC.CountBits();
+ }
+ else {
+ AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+ const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
+ allclusters = clusterTPC.CountBits();
+ }
+ clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
+ }
return clusterRatio;
-}
-//______________________________________________________
+}
+//___________________________________________________
void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
//
// Get impact parameter
//
+
+ const Double_t kBeampiperadius=3.;
+ Double_t dcaD[2]={-999.,-999.},
+ covD[3]={-999.,-999.,-999.};
TClass *type = track->IsA();
if(type == AliESDtrack::Class()){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
esdtrack->GetImpactParameters(radial, z);
}
else if(type == AliAODTrack::Class()){
- AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
- if(aodtrack){
- Double_t xyz[3];
- aodtrack->XYZAtDCA(xyz);
- //printf("xyz[0] %f, xyz[1] %f, xyz[2] %f\n",xyz[0], xyz[1],xyz[2]);
- if(xyz[0]==-999. || xyz[1]==-999. || xyz[2]==-999.){
- if(!fEvent) {
- z = xyz[2];
- radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
- //printf("No event\n");
- }
- else {
- // Propagate the global track to the DCA.
- const AliVVertex *vAOD = fEvent->GetPrimaryVertex();
- Double_t PosAtDCA[2] = {-999,-999};
- Double_t covar[3] = {-999,-999,-999};
- AliAODTrack *copyaodtrack = new AliAODTrack(*aodtrack);
- if(copyaodtrack->PropagateToDCA(vAOD,fEvent->GetMagneticField(),100.,PosAtDCA,covar)) {
- z = PosAtDCA[1];
- radial = PosAtDCA[0];
- //printf("Propagate z %f and radial %f\n",z,radial);
- }
- if(copyaodtrack) delete copyaodtrack;
- }
- }
- else {
- z = xyz[2];
- radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
- }
+
+ //case of AOD tracks
+ AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+ if(!aodevent) {
+ AliDebug(1, "No aod event available\n");
+ return;
+ }
+
+ //Case ESD track: take copy constructor
+ AliAODTrack *aodtrack = NULL;
+ AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+ if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+ AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
+ if(!vtxAODSkip) return;
+ AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+ if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+ radial = dcaD[0];
+ z = dcaD[1];
}
+ //if(vtxAODSkip) delete vtxAODSkip;
+ if(aodtrack) delete aodtrack;
}
}
//______________________________________________________
if(type == AliESDtrack::Class()){
AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
if(!esdtrack) return kFALSE;
- if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
- else return kTRUE;
+ if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
+ else return kFALSE;
}
else if(type == AliAODTrack::Class()){
return kFALSE;
}
+
+//______________________________________________________
+Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
+ //
+ // Get TRDchi2
+ //
+ Int_t ntls(0);
+ TClass *type = track->IsA();
+ if(type == AliESDtrack::Class()){
+ AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+ ntls = esdtrack->GetTRDntracklets();
+ return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
+ }
+ else if(type == AliAODTrack::Class()){
+ AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+ if(aodtrack){
+ return 999.;
+ }
+ }
+
+ return 999.;
+
+}
+
//______________________________________________________
Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
//
return 0;
}
//______________________________________________________
-void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
- //
- // Get HFE impact parameter (with recalculated primary vertex)
- //
- dcaxy=0;
- dcansigmaxy=0;
+void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
+ //
+ // Get HFE impact parameter (with recalculated primary vertex)
+ //
+ dcaxy=0;
+ dcansigmaxy=0;
if(!fEvent){
AliDebug(1, "No Input event available\n");
return;
}
- const Double_t kBeampiperadius=3.;
TString type = track->IsA()->GetName();
- Double_t dcaD[2]={-999.,-999.};
- Double_t covD[3]={-999.,-999.,-999.};
- Float_t dcaF[2]={-999.,-999.};
- Float_t covF[3]={-999.,-999.,-999.};
-
- AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
- if(!esdevent) {
- AliDebug(1, "No esd event available\n");
- return;
+ const Double_t kBeampiperadius=3.;
+ Double_t dcaD[2]={-999.,-999.},
+ covD[3]={-999.,-999.,-999.};
+ Bool_t isRecalcVertex(kFALSE);
+
+ if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+ //case of ESD tracks
+ AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+ if(!esdevent) {
+ AliDebug(1, "No esd event available\n");
+ return;
+ }
+
+ const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+ if(!vtxESDSkip) return;
+
+ //case ESD track: take copy constructor
+ const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
+ if(tmptrack){
+
+ if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+ vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+ isRecalcVertex = kTRUE;
+ }
+
+ if(vtxESDSkip){
+ AliESDtrack esdtrack(*tmptrack);
+ fMagField = fEvent->GetMagneticField();
+ if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
+ dcaxy = dcaD[0];
+ if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+ }
+ if(isRecalcVertex) delete vtxESDSkip;
+ }
+ }
}
- const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
- if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
- // recalculate primary vertex for peri. and pp
- AliVertexerTracks vertexer(fEvent->GetMagneticField());
- vertexer.SetITSMode();
- vertexer.SetMinClusters(4);
- Int_t skipped[2];
- skipped[0] = track->GetID();
- vertexer.SetSkipTracks(1,skipped);
-
- //----diamond constraint-----------------------------
- vertexer.SetConstraintOn();
- Float_t diamondcovxy[3];
- esdevent->GetDiamondCovXY(diamondcovxy);
- Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
- Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
- AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
- vertexer.SetVtxStart(diamond);
- delete diamond; diamond=NULL;
- //----------------------------------------------------
-
- vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
-
- // Getting the DCA
- // Propagation always done on a working copy to not disturb the track params of the original track
- AliESDtrack *esdtrack = NULL;
- if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
- // Case ESD track: take copy constructor
- AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
- if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
- } else {
- // Case AOD track: take different constructor
- esdtrack = new AliESDtrack(track);
- }
- if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
- // protection
- dcaxy = dcaD[0];
- if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
- if(!covD[0]) dcansigmaxy = -49.;
- }
- delete esdtrack;
- delete vtxESDSkip;
- }
- else{
- AliESDtrack *esdtrack = NULL;
- if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
- // Case ESD track: take copy constructor
- AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
- if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
- } else {
- // Case AOD track: take different constructor
- esdtrack = new AliESDtrack(track);
- }
- if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
- dcaxy = dcaF[0];
- if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
- if(!covF[0]) dcansigmaxy = -49.;
- delete esdtrack;
+ else {
+ //case of AOD tracks
+ AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+ if(!aodevent) {
+ AliDebug(1, "No aod event available\n");
+ return;
+ }
+
+ AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
+ if(!vtxAODSkip) return;
+
+ //Case ESD track: take copy constructor
+ const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
+ if(tmptrack){
+
+ if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+ vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+ isRecalcVertex = kTRUE;
+ }
+ if(vtxAODSkip){
+ AliAODTrack aodtrack(*tmptrack);
+ AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
+ fMagField = aodevent->GetMagneticField();
+ if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
+ dcaxy = dcaD[0];
+ if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+ }
+ if(isRecalcVertex) delete vtxAODSkip;
+ }
+ }
}
}
//______________________________________________________
-void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
+void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
//
// Get HFE impact parameter (with recalculated primary vertex)
//
}
const Double_t kBeampiperadius=3.;
TString type = track->IsA()->GetName();
- Float_t dcaF[2]={-999.,-999.};
- Float_t covF[3]={-999.,-999.,-999.};
+ Bool_t isRecalcVertex(kFALSE);
- AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
- if(!esdevent) {
- AliDebug(1, "No esd event available\n");
- return;
+ if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+ //case of ESD tracks
+ AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+ if(!esdevent) {
+ AliDebug(1, "No esd event available\n");
+ return;
+ }
+
+ // Check whether primary vertex is available
+ const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+ if(!vtxESDSkip) return;
+
+ const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
+ if(tmptrack){
+ //case ESD track: take copy constructor
+
+ if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+ vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+ isRecalcVertex = kTRUE;
+ }
+ if(vtxESDSkip){
+ AliESDtrack esdtrack(*tmptrack);
+ fMagField = fEvent->GetMagneticField();
+ esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
+
+ if(isRecalcVertex) delete vtxESDSkip;
+ }
+ }
}
- const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
- if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
- // recalculate primary vertex for peri. and pp
- AliVertexerTracks vertexer(fEvent->GetMagneticField());
- vertexer.SetITSMode();
- vertexer.SetMinClusters(4);
- Int_t skipped[2];
- skipped[0] = track->GetID();
- vertexer.SetSkipTracks(1,skipped);
-
- //----diamond constraint-----------------------------
- vertexer.SetConstraintOn();
- Float_t diamondcovxy[3];
- esdevent->GetDiamondCovXY(diamondcovxy);
- Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
- Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
- AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
- vertexer.SetVtxStart(diamond);
- delete diamond; diamond=NULL;
- //----------------------------------------------------
-
- vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
-
- // Getting the DCA
- // Propagation always done on a working copy to not disturb the track params of the original track
- AliESDtrack *esdtrack = NULL;
- if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
- // Case ESD track: take copy constructor
- AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
- if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
- } else {
- // Case AOD track: take different constructor
- esdtrack = new AliESDtrack(track);
- }
- if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
- delete esdtrack;
- delete vtxESDSkip;
- }
- else{
- AliESDtrack *esdtrack = NULL;
- if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
- // Case ESD track: take copy constructor
- AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
- if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
- } else {
- // Case AOD track: take different constructor
- esdtrack = new AliESDtrack(track);
- }
- if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
- dcaD[0] = dcaF[0];
- dcaD[1] = dcaF[1];
- covD[0] = covF[0];
- covD[1] = covF[1];
- covD[2] = covF[2];
- delete esdtrack;
+ else {
+ //case of AOD tracks
+ AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+ if(!aodevent) {
+ AliDebug(1, "No aod event available\n");
+ return;
+ }
+
+ // Check whether primary vertex is available
+ AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
+ if(!vtxAODSkip) return;
+
+ //Case ESD track: take copy constructor
+ const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
+ if(tmptrack){
+
+ if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+ vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+ isRecalcVertex = kTRUE;
+ }
+ if(vtxAODSkip){
+ AliAODTrack aodtrack(*tmptrack);
+ AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
+ fMagField = aodevent->GetMagneticField();
+ etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
+
+ if(isRecalcVertex) delete vtxAODSkip;
+ }
+ }
}
}
//______________________________________________________
-void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
+void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
//
// Get HFE impact parameter cut(pt dependent)
//
- TString type = track->IsA()->GetName();
- if(!type.CompareTo("AliESDtrack")){
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
- if(!esdtrack) return;
- Double_t pt = esdtrack->Pt();
- hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
- hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
- }
+ Double_t pt = track->Pt();
+ hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
+ hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
}
//______________________________________________________
-void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
+void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
//
// Get max impact parameter cut r (pt dependent)
//
+ Double_t pt = track->Pt();
+ if(pt > 0.15) {
+ maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
+ }
+ else maximpactRcut = 9999999999.0;
+}
+//______________________________________________________
+void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
+ //
+ // TOF matching
+ //
+
TString type = track->IsA()->GetName();
if(!type.CompareTo("AliESDtrack")){
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+ const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
if(!esdtrack) return;
- Double_t pt = esdtrack->Pt();
- if(pt > 0.15) {
- maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
+ tofsignalDx = esdtrack->GetTOFsignalDx();
+ tofsignalDz = esdtrack->GetTOFsignalDz();
+ }
+
+}
+
+//______________________________________________________
+Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
+ //
+ // Check if every ITS layer, which has a module which is alive, also
+ // has an ITS cluster
+ //
+ Bool_t patternOK(kTRUE);
+ Int_t status(0);
+ for(Int_t ily = 0; ily < 6; ily++){
+ status = GetITSstatus(track, ily);
+ if(CheckITSstatus(status)){
+ // pixel alive, check whether layer has a cluster
+ if(!TESTBIT(track->GetITSClusterMap(),ily)){
+ // No cluster even though pixel is alive - reject track
+ patternOK = kFALSE;
+ break;
+ }
}
- else maximpactRcut = 9999999999.0;
}
+ return patternOK;
+}
+
+//---------------------------------------------------------------------------
+const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
+ //
+ // This method returns a primary vertex without the daughter tracks of the
+ // candidate and it recalculates the impact parameters and errors for ESD tracks.
+ //
+ // The output vertex is created with "new".
+ //
+
+ //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+
+ AliVertexerTracks vertexer(fEvent->GetMagneticField());
+ vertexer.SetITSMode();
+ vertexer.SetMinClusters(4);
+ Int_t skipped[2];
+ skipped[0] = track->GetID();
+ vertexer.SetSkipTracks(1,skipped);
+
+ //diamond constraint
+ vertexer.SetConstraintOn();
+ Float_t diamondcovxy[3];
+ esdevent->GetDiamondCovXY(diamondcovxy);
+ Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
+ Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+ AliESDVertex diamond(pos,cov,1.,1);
+ vertexer.SetVtxStart(&diamond);
+
+ const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
+
+ return vtxESDSkip;
+}
+
+//---------------------------------------------------------------------------
+AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
+ //
+ // This method returns a primary vertex without the daughter tracks of the
+ // candidate and it recalculates the impact parameters and errors for AOD tracks.
+ // The output vertex is created with "new".
+ //
+
+ AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+ if(!vtxAOD) return 0;
+ TString title=vtxAOD->GetTitle();
+ if(!title.Contains("VertexerTracks")) return 0;
+
+ AliVertexerTracks vertexer(aod->GetMagneticField());
+
+ vertexer.SetITSMode();
+ vertexer.SetMinClusters(3);
+ vertexer.SetConstraintOff();
+
+ if(title.Contains("WithConstraint")) {
+ Float_t diamondcovxy[3];
+ aod->GetDiamondCovXY(diamondcovxy);
+ Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
+ Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+ AliESDVertex diamond(pos,cov,1.,1);
+ vertexer.SetVtxStart(&diamond);
+ }
+ Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
+ Int_t id = (Int_t)track->GetID();
+ if(!(id<0)) skipped[0] = id;
+
+ /*// exclude tracks with global constrained parameters
+ Int_t nTracks=aod->GetNumberOfTracks();
+ for(Int_t i=0; i<nTracks; i++){
+ t = aod->GetTrack(i);
+ if(t->TestFilterMask(512)){
+ id = (Int_t)t->GetID();
+ skipped[nTrksToSkip++] = id;
+ }
+ }*/
+
+ vertexer.SetSkipTracks(1,skipped);
+ AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
+
+ if(!vtxESDNew) return 0;
+ if(vtxESDNew->GetNContributors()<=0) {
+ delete vtxESDNew; vtxESDNew=NULL;
+ return 0;
+ }
+
+ // convert to AliAODVertex
+ Double_t pos[3],cov[6],chi2perNDF;
+ vtxESDNew->GetXYZ(pos); // position
+ vtxESDNew->GetCovMatrix(cov); //covariance matrix
+ chi2perNDF = vtxESDNew->GetChi2toNDF();
+ delete vtxESDNew; vtxESDNew=NULL;
+
+ AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
+
+ return vtxAODNew;
}