#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMCParticle.h"
+#include "AliStack.h"
#include "AliPIDResponse.h"
+#include "AliTrackReference.h"
#include "AliVEvent.h"
+#include "AliHFEpidTPC.h"
+#include "AliHFEpidTRD.h"
+#include "AliHFEmcQA.h"
#include "TTreeStream.h"
#include "AliHFEdebugTreeTask.h"
AliAnalysisTaskSE(),
fTrackCuts(NULL),
fSignalCuts(NULL),
+ fTRDpid(NULL),
+ fTPCpid(NULL),
+ fExtraCuts(NULL),
fNclustersTPC(70),
fNclustersTPCPID(0),
fNclustersITS(2),
fFilename("HFEtree.root"),
- fDebugTree(NULL)
+ fDebugTree(NULL),
+ fNparents(-1)
{
}
AliAnalysisTaskSE(name),
fTrackCuts(NULL),
fSignalCuts(NULL),
+ fTRDpid(NULL),
+ fTPCpid(NULL),
+ fExtraCuts(NULL),
fNclustersTPC(70),
fNclustersTPCPID(0),
fNclustersITS(2),
fFilename("HFEtree.root"),
- fDebugTree(NULL)
+ fDebugTree(NULL),
+ fNparents(-1)
{
-
+ fTRDpid = new AliHFEpidTRD("QAtrdPID");
+ fTPCpid = new AliHFEpidTPC("QAtpcPID");
}
AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
- if(fDebugTree) delete fDebugTree;
+ if(fDebugTree) delete fDebugTree;
+ if(fTRDpid) delete fTRDpid;
+ if(fTPCpid) delete fTPCpid;
}
void AliHFEdebugTreeTask::UserCreateOutputObjects(){
fTrackCuts->SetUseMixedVertex(kTRUE);
fTrackCuts->SetVertexRange(10.);
fTrackCuts->Initialize();
+
+ fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+
}
void AliHFEdebugTreeTask::UserExec(Option_t *){
AliError("No Handler");
}
if(!pid){
- AliError("No PID response\n");
+ AliError("No PID response");
return;
}
+ if(!fInputEvent) {
+ AliError("No Input event");
+ return;
+ }
+
+ AliESDtrack copyTrack;
fTrackCuts->SetRecEvent(fInputEvent);
// Cut event
if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
- AliDebug(1, "Event flagged as pileup\n");
+ AliDebug(1, "Event flagged as pileup\n");
return;
}
if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
fSignalCuts->SetMCEvent(fMCEvent);
}
+ // MC get stack
+ AliStack* stack = 0x0;
+ if(mcthere){
+ stack = fMCEvent->Stack();
+ if(!stack) AliError("No Stack");
+ }
+
// Get Primary Vertex
const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
+ Double_t vtx[3];
+ vertex->GetXYZ(vtx);
+ Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
// Get centrality
Float_t centrality = -1.;
- TString beamtype = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetESDRun()->GetBeamType();
+ AliESDEvent *event = (dynamic_cast<AliESDEvent *>(fInputEvent));
+ if(!event) return;
+ TString beamtype = event->GetBeamType();
+ //printf("Beam type %s\n",(const char*)beamtype);
if(!beamtype.CompareTo("Pb-Pb") || !beamtype.CompareTo("A-A")){
// Heavy ion run
AliDebug(1, "Heavy-Ion event\n");
centrality = hicent->GetCentralityPercentile("V0M");
}
+ if(!fExtraCuts){
+ fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+ }
+ fExtraCuts->SetRecEventInfo(event);
+
+ // Store event selection variables
+ (*fDebugTree) << "EventDebug"
+ << "Centrality=" << centrality
+ << "VertexZ=" << vtx[2]
+ << "NumberOfContributors=" << ncontrib
+ << "run=" << run
+ << "\n";
+
// Common variables
Double_t charge, eta, phi, momentum, transversemomentum;
Int_t source;
// Get Production Vertex in radial direction
Double_t vx = mcpart->Particle()->Vx(),
- vy = mcpart->Particle()->Vy();
+ vy = mcpart->Particle()->Vy();
Double_t productionVertex = TMath::Sqrt(vx*vx+vy*vy);
// Get Mother PDG code of the particle
else if(TMath::Abs(pdg) == 11) source = 4;
else source = 5;
- (*fDebugTree) << "MCDebug"
+ // Momemntum at the inner wall of the TPC
+ Double_t pTPC = 0., ptTPC = 0.;
+ AliTrackReference *ref = FindTrackReference(mcpart, 80, 270, AliTrackReference::kTPC);
+ if(ref){
+ pTPC = ref->P();
+ ptTPC = ref->Pt();
+ }
+
+
+
+ (*fDebugTree) << "MCDebug"
<< "centrality=" << centrality
<< "MBtrigger=" << isMBTrigger
<< "CentralTrigger=" << isCentralTrigger
}
AliESDtrack *track;
+ Double_t mcp, mcpt, mcptTPC, mcpTPC; // MC Variables added to the debug tree
for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
// fill the tree
track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
// Cut track (Only basic track cuts)
if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
// Debug streaming of PID-related quantities
+ copyTrack.~AliESDtrack();
+ new(©Track) AliESDtrack(*track);
+ if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(©Track, AliHFEpidObject::kESDanalysis); // Apply Eta Correction on copy track
Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
- Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
- if(TMath::Abs(nSigmaTOF) > 5) continue;
+ Double_t nSigmaTPC = pid->NumberOfSigmasTPC(©Track, AliPID::kElectron);
+ //if(TMath::Abs(nSigmaTOF) > 5) continue;
// we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
- Double_t tPCdEdx = track->GetTPCsignal();
+ Double_t tPCdEdx = copyTrack.GetTPCsignal();
// Signal, source and MCPID
Bool_t signal = kTRUE;
source = 5;
- if(mcthere){
+ mcp = mcpt = mcpTPC = mcptTPC = 0.;
+
+
+
+ Double_t bgcategory = 0.;
+ Int_t mArr = -1;
+ Int_t mesonID = -999;
+ Double_t xr[3]={-999,-999,-999};
+ Double_t eR=-999;
+ Double_t eZ=-999;
+ Double_t unique=-999;
+ Double_t mesonunique=-999;
+ Double_t mesonR=-999;
+ Double_t mesonZ=-999;
+ Double_t mesonMomPdg=-999;
+ Double_t mesonMomPt=-999;
+ Double_t mesonGMomPdg=-999;
+ Double_t mesonGGMomPdg=-999;
+ Double_t mesonPt = -999;
+ Double_t mceta = -999;
+ Double_t mcphi = -999;
+ Int_t mcpdg;
+
+ if(mcthere){
// Signal
AliMCParticle *mctrack;
if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))){
else if(fSignalCuts->IsNonHFElectron(track)) source = 3;
else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)) source = 4;
else source = 5;
+
+ if(!mctrack) continue;
+
+ // Kinematics
+ mcpt = mctrack->Pt();
+ mcp = mctrack->P();
+ mceta = mctrack->Eta();
+ mcphi = mctrack->Phi();
+ mcpdg = mctrack->Particle()->GetPdgCode();
+
+ AliTrackReference *ref = FindTrackReference(mctrack, 80, 270, AliTrackReference::kTPC);
+ if(ref){
+ mcpTPC = ref->P();
+ mcptTPC = ref->Pt();
+ }
+
+
+ TParticle *mctrack1 = mctrack->Particle();
+ mesonID=GetElecSourceMC(mctrack1);
+ if(mesonID==AliHFEmcQA::kGammaPi0 || mesonID==AliHFEmcQA::kPi0) mArr=0; //pion
+ else if(mesonID==AliHFEmcQA::kGammaEta || mesonID==AliHFEmcQA::kEta) mArr=1; //eta
+ else if(mesonID==AliHFEmcQA::kGammaOmega || mesonID==AliHFEmcQA::kOmega) mArr=2; //omega
+ else if(mesonID==AliHFEmcQA::kGammaPhi || mesonID==AliHFEmcQA::kPhi) mArr=3; //phi
+ else if(mesonID==AliHFEmcQA::kGammaEtaPrime || mesonID==AliHFEmcQA::kEtaPrime) mArr=4; //etaprime
+ else if(mesonID==AliHFEmcQA::kGammaRho0 || mesonID==AliHFEmcQA::kRho0) mArr=5; //rho
+
+ mctrack->XvYvZv(xr);
+
+ eR= TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+ eZ = xr[2];
+ TParticle *mctrackt = mctrack->Particle();
+ unique=mctrackt->GetUniqueID();
+
+ AliMCParticle *mctrackmother = NULL;
+ AliMCParticle *mctrackmother2 = NULL;
+
+ if(!(mArr<0)){
+ if(mesonID>=AliHFEmcQA::kGammaPi0) { // conversion electron, be careful with the enum odering
+ Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonPt = mctrackmother->Pt(); //meson pt
+ bgcategory = 1.;
+ mctrackmother->XvYvZv(xr);
+ mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+ mesonZ = xr[2];
+
+ mctrackt = mctrackmother->Particle();
+ if(mctrackt){
+ mesonunique = mctrackt->GetUniqueID();
+ }
+
+ Int_t glabel2=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+ if((mctrackmother2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel2)))){
+ mesonMomPdg=mctrackmother2->PdgCode();
+ mesonMomPt=mctrackmother2->Pt();
+ }
+
+ if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+ bgcategory = 2.;
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonMomPdg=mctrackmother->PdgCode();
+ mesonMomPt=mctrackmother->Pt();
+ if(TMath::Abs(mctrackmother->PdgCode())==310){
+ bgcategory = 3.;
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonGMomPdg=mctrackmother->PdgCode();
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonGGMomPdg=mctrackmother->PdgCode();
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ else{ // nonHFE except for the conversion electron
+ Int_t glabel=TMath::Abs(mctrack->GetMother());
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonPt = mctrackmother->Pt(); //meson pt
+ bgcategory = -1.;
+ mctrackmother->XvYvZv(xr);
+ mesonR = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+ mesonZ = xr[2];
+
+ mctrackt = mctrackmother->Particle();
+ if(mctrackt){
+ mesonunique = mctrackt->GetUniqueID();
+ }
+ if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+ bgcategory = -2.;
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonMomPdg=mctrackmother->PdgCode();
+ mesonMomPt=mctrackmother->Pt();
+ if(TMath::Abs(mctrackmother->PdgCode())==310){
+ bgcategory = -3.;
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonGMomPdg=mctrackmother->PdgCode();
+ glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+ if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+ mesonGGMomPdg=mctrackmother->PdgCode();
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+
}
// Get V0 tag (if available)
Int_t v0pid = -1;
phi = track->Phi();
momentum = track->P() * charge;
transversemomentum = track->Pt() * charge;
+ Double_t momentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : 0.;
+ Double_t transversemomentumTPC = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->Pt() : 0.;
// ITS number of clusters
UChar_t nclustersITS = track->GetITSclusters(NULL);
Double_t chi2matching = track->GetChi2TPCConstrainedVsGlobal(dynamic_cast<const AliESDVertex *>(vertex));
Double_t chi2PerClusterITS = 0.0;
if (nclustersITS != 0) chi2PerClusterITS = track->GetITSchi2() / Float_t(nclustersITS);
// TPC number of clusters (different definitions)
- UChar_t nclustersTPC = track->GetTPCncls();
+ UChar_t nclustersTPCfit = track->GetTPCNcls();
+ UChar_t nclustersTPCall = 0;
+ const TBits &clusterTPC = track->GetTPCClusterMap();
+ nclustersTPCall = clusterTPC.CountBits();
UChar_t nclustersTPCPID = track->GetTPCsignalN();
UChar_t nfindableTPC = track->GetTPCNclsF();
- Double_t clusterRatioTPC = 0.0;
- if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPC = static_cast<Double_t>(nclustersTPC)/static_cast<Double_t>(nfindableTPC);
+ Double_t clusterRatioTPCfit = 0.0;
+ if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC);
+ Double_t clusterRatioTPCall = 0.0;
+ if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC);
UChar_t nclustersTPCshared = 0;
Float_t ncrossedRowsTPC = track->GetTPCCrossedRows();
const TBits &sharedTPC = track->GetTPCSharedMap();
for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
// TRD clusters and tracklets
UChar_t nclustersTRD = track->GetTRDncls();
- UChar_t ntracklets = track->GetTRDntrackletsPID();
+ UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
// ITS and TRD acceptance maps
- UChar_t hasClusterITS[6], hasTrackletTRD[6];
+ UChar_t hasClusterITS[6], statusITS[6], hasTrackletTRD[6];
UChar_t itsPixel = track->GetITSClusterMap();
- for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
+ for(Int_t icl = 0; icl < 6; icl++){
+ hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
+ if(CheckITSstatus(track, icl)) statusITS[icl] = 1;
+ else statusITS[icl] = 0;
+ }
+ Double_t trddEdxSum[6];
+ for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
for(Int_t itl = 0; itl < 6; itl++){
Int_t nSliceNonZero = 0;
+ trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
for(Int_t islice = 0; islice < 8; islice++){
if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
}
track->GetTRDpid(pidprobs);
Double_t likeEleTRD = pidprobs[0];
Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
+ Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
+ Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);
+
// DCA
Float_t b[2] = {0.,0.};
Float_t bCov[3] = {0.,0.,0.};
if(bCov[0]>0) dcaSR = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
if(bCov[2]>0) dcaSZ = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
Double_t dcaS = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
- // Vertex
- Double_t vtx[3];
- vertex->GetXYZ(vtx);
- Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
+
+ // HFE DCA
+ Double_t hfeb[2] = {-99.,-99.};
+ Double_t hfebCov[3] = {-999.,-999.,-999.};
+ fExtraCuts->GetHFEImpactParameters(track, hfeb, hfebCov);
+
+ Double_t tofdx= -999.0;
+ Double_t tofdz= -999.0;
+ tofdx=track->GetTOFsignalDx();
+ tofdz=track->GetTOFsignalDz();
+
+ // TOF track status
+ UInt_t status = 0;
+ status = track->GetStatus();
+ Bool_t hasTOFout = status&AliESDtrack::kTOFout;
+ Bool_t hasTOFtime = status&AliESDtrack::kTIME;
+ Bool_t hasTOFpid = status&AliESDtrack::kTOFpid;
+ Bool_t hasgoodTOF = kFALSE;
+ if (hasTOFout && hasTOFtime && hasTOFpid) hasgoodTOF = kTRUE;
+
+ // TRD track status
+ Bool_t hasTRDin = status&AliESDtrack::kTRDin;
+
+
+ // TOF mismatch (particle spectra group)
+ Int_t mismatchlevel=0; // default value; in data always 0
+ if(mcthere){
+ Int_t tofLabel[3];
+ track->GetTOFLabel(tofLabel);
+ if(TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) mismatchlevel=1;
+ TParticle *matchedTrack = stack->Particle(TMath::Abs(tofLabel[0]));
+ if(TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
+ {
+ if(mismatchlevel==1) mismatchlevel=3;
+ else mismatchlevel=2;
+ }
+ }
+
+
// Fill Tree
(*fDebugTree) << "PIDdebug"
<< "centrality=" << centrality
<< "v0pid=" << v0pid
<< "run=" << run
<< "p=" << momentum
+ << "ptpc=" << momentumTPC
<< "pt=" << transversemomentum
+ << "pttpc=" << transversemomentumTPC
+ << "mcp=" << mcp
+ << "mcpt=" << mcpt
+ << "mcpTPC=" << mcpTPC
+ << "mcptTPC=" << mcptTPC
+ << "mceta=" << mceta
+ << "mcphi=" << mcphi
+ << "mcpdg=" << mcpdg
<< "eta=" << eta
<< "phi=" << phi
- << "ntracklets=" << ntracklets
- << "nclustersTPC=" << nclustersTPC
+ << "ntracklets=" << ntrackletsTRDPID
+ << "nclustersTPC=" << nclustersTPCfit
+ << "nclustersTPCall=" << nclustersTPCall
<< "nclustersTPCPID=" << nclustersTPCPID
<< "nclustersTPCshared=" << nclustersTPCshared
<< "ncrossedRowsTPC=" << ncrossedRowsTPC
- << "clusterRatioTPC=" << clusterRatioTPC
+ << "clusterRatioTPC=" << clusterRatioTPCfit
+ << "clusterRatioTPCall=" << clusterRatioTPCall
<< "nclustersITS=" << nclustersITS
<< "nclusters=" << nclustersTRD
<< "chi2matching=" << chi2matching
- << "chi2PerClusterITS=" << chi2PerClusterITS
+ << "chi2PerClusterITS=" << chi2PerClusterITS
<< "its0=" << hasClusterITS[0]
<< "its1=" << hasClusterITS[1]
<< "its2=" << hasClusterITS[2]
<< "its3=" << hasClusterITS[3]
<< "its4=" << hasClusterITS[4]
<< "its5=" << hasClusterITS[5]
+ << "statusITS0=" << statusITS[0]
+ << "statusITS1=" << statusITS[1]
+ << "statusITS2=" << statusITS[2]
+ << "statusITS3=" << statusITS[3]
+ << "statusITS4=" << statusITS[4]
+ << "statusITS5=" << statusITS[5]
<< "trd0=" << hasTrackletTRD[0]
<< "trd1=" << hasTrackletTRD[1]
<< "trd2=" << hasTrackletTRD[2]
<< "trd3=" << hasTrackletTRD[3]
<< "trd4=" << hasTrackletTRD[4]
<< "trd5=" << hasTrackletTRD[5]
+ << "TRDdEdxl0=" << trddEdxSum[0]
+ << "TRDdEdxl1=" << trddEdxSum[1]
+ << "TRDdEdxl2=" << trddEdxSum[2]
+ << "TRDdEdxl3=" << trddEdxSum[3]
+ << "TRDdEdxl4=" << trddEdxSum[4]
+ << "TRDdEdxl5=" << trddEdxSum[5]
<< "TOFsigmaEl=" << nSigmaTOF
<< "TPCsigmaEl=" << nSigmaTPC
<< "TPCdEdx=" << tPCdEdx
<< "TRDlikeEl=" << likeEleTRD
<< "TRDlikeEln=" << likeEleTRDn
+ << "trdtruncmean1=" << trdtruncmean1
+ << "trdtruncmean2=" << trdtruncmean2
<< "dcaR=" << b[0]
<< "dcaZ=" << b[1]
<< "dca=" << dca
<< "dcaSR=" << dcaSR
<< "dcaSZ=" << dcaSZ
<< "dcaS=" << dcaS
+ << "hfedcaR=" << hfeb[0]
+ << "hfedcaZ=" << hfeb[1]
+ << "hfedcacovR=" << hfebCov[0]
+ << "hfedcacovZ=" << hfebCov[2]
<< "vx=" << vtx[0]
<< "vy=" << vtx[1]
- << "vz=" << vtx[2]
+ << "vz=" << vtx[2]
+ << "tofdx=" << tofdx
+ << "tofdz=" << tofdz
+ << "statusTOFtracking=" << hasgoodTOF
+ << "TOFmismatchlevel=" << mismatchlevel
+ << "statusTRDtracking=" << hasTRDin
<< "ncontrib=" << ncontrib
+ << "mesonID=" << mesonID
+ << "eR=" << eR
+ << "mesonR=" << mesonR
+ << "eZ=" << eZ
+ << "mesonZ=" << mesonZ
+ << "unique=" << unique
+ << "mesonunique=" << mesonunique
+ << "bgcategory=" << bgcategory
+ << "mesonpt=" << mesonPt
+ << "mesonMomPdg=" << mesonMomPdg
+ << "mesonGMomPdg=" << mesonGMomPdg
+ << "mesonGGMomPdg=" << mesonGGMomPdg
+ << "mesonMomPt=" << mesonMomPt
<< "\n";
}
}
void AliHFEdebugTreeTask::SetFileName(const char *filename){ fFilename = filename; }
+//___________________________________________________________
+AliTrackReference *AliHFEdebugTreeTask::FindTrackReference(AliMCParticle *track, Float_t minRadius, Float_t maxRadius, Int_t detectorID)
+{
+ //
+ // Find the track reference
+ //
+ AliTrackReference *ref = NULL, *reftmp;
+ Float_t radius;
+ for(Int_t iref = 0; iref < track->GetNumberOfTrackReferences(); iref++){
+ reftmp = track->GetTrackReference(iref);
+ if(reftmp->DetectorId() != detectorID) continue;
+ radius = reftmp->R();
+ if(radius >= minRadius && radius < maxRadius){
+ ref = reftmp;
+ break;
+ }
+ if(radius > maxRadius) break;
+ }
+ return ref;
+}
+
+//__________________________________________
+Int_t AliHFEdebugTreeTask::GetElecSourceMC(TParticle * const mcpart)
+{
+ // decay particle's origin
+
+ if(!mcpart){
+ AliDebug(1, "no mcparticle, return\n");
+ return -1;
+ }
+
+ if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return AliHFEmcQA::kMisID;
+
+ Int_t origin = -1;
+ Bool_t isFinalOpenCharm = kFALSE;
+
+ Int_t iLabel = mcpart->GetFirstMother();
+ if (iLabel<0){
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ AliMCParticle *mctrack = NULL;
+ Int_t tmpMomLabel=0;
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
+ TParticle *partMother = mctrack->Particle();
+ TParticle *partMotherCopy = mctrack->Particle();
+ Int_t maPdgcode = partMother->GetPdgCode();
+
+ // if the mother is charmed hadron
+ if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kCharm || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kCharm ) {
+
+ for (Int_t i=0; i<fNparents; i++){
+ if (abs(maPdgcode)==fParentSelect[0][i]){
+ isFinalOpenCharm = kTRUE;
+ }
+ }
+ if (!isFinalOpenCharm) return -1;
+
+ // iterate until you find B hadron as a mother or become top ancester
+ for (Int_t i=1; i<fgkMaxIter; i++){
+
+ Int_t jLabel = partMother->GetFirstMother();
+ if (jLabel == -1){
+ origin = AliHFEmcQA::kDirectCharm;
+ return origin;
+ }
+ if (jLabel < 0){ // safety protection
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ // if there is an ancester
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
+ TParticle* grandMa = mctrack->Particle();
+ Int_t grandMaPDG = grandMa->GetPdgCode();
+
+ for (Int_t j=0; j<fNparents; j++){
+ if (abs(grandMaPDG)==fParentSelect[1][j]){
+ origin = AliHFEmcQA::kBeautyCharm;
+ return origin;
+ }
+ }
+
+ partMother = grandMa;
+ } // end of iteration
+ } // end of if
+ else if ( (int(abs(maPdgcode)/100.)%10) == AliHFEmcQA::kBeauty || (int(abs(maPdgcode)/1000.)%10) == AliHFEmcQA::kBeauty ) {
+ for (Int_t i=0; i<fNparents; i++){
+ if (abs(maPdgcode)==fParentSelect[1][i]){
+ origin = AliHFEmcQA::kDirectBeauty;
+ return origin;
+ }
+ }
+ } // end of if
+ else if ( abs(maPdgcode) == 22 ) { //conversion
+
+ tmpMomLabel = partMotherCopy->GetFirstMother();
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
+ partMother = mctrack->Particle();
+ maPdgcode = partMother->GetPdgCode();
+ if ( abs(maPdgcode) == 111 ) {
+ origin = AliHFEmcQA::kGammaPi0;
+ return origin;
+ }
+ else if ( abs(maPdgcode) == 221 ) {
+ origin = AliHFEmcQA::kGammaEta;
+ return origin;
+ }
+ else if ( abs(maPdgcode) == 223 ) {
+ origin = AliHFEmcQA::kGammaOmega;
+ return origin;
+ }
+ else if ( abs(maPdgcode) == 333 ) {
+ origin = AliHFEmcQA::kGammaPhi;
+ return origin;
+ }
+ else if ( abs(maPdgcode) == 331 ) {
+ origin = AliHFEmcQA::kGammaEtaPrime;
+ return origin;
+ }
+ else if ( abs(maPdgcode) == 113 ) {
+ origin = AliHFEmcQA::kGammaRho0;
+ return origin;
+ }
+ else origin = AliHFEmcQA::kElse;
+ //origin = kGamma; // finer category above
+ return origin;
+
+ } // end of if
+ else if ( abs(maPdgcode) == 111 ) {
+ origin = AliHFEmcQA::kPi0;
+ return origin;
+ } // end of if
+ else if ( abs(maPdgcode) == 221 ) {
+ origin = AliHFEmcQA::kEta;
+ return origin;
+ } // end of if
+ else if ( abs(maPdgcode) == 223 ) {
+ origin = AliHFEmcQA::kOmega;
+ return origin;
+ } // end of if
+ else if ( abs(maPdgcode) == 333 ) {
+ origin = AliHFEmcQA::kPhi;
+ return origin;
+ } // end of if
+ else if ( abs(maPdgcode) == 331 ) {
+ origin = AliHFEmcQA::kEtaPrime;
+ return origin;
+ } // end of if
+ else if ( abs(maPdgcode) == 113 ) {
+ origin = AliHFEmcQA::kRho0;
+ return origin;
+ } // end of if
+ else{
+ origin = AliHFEmcQA::kElse;
+ }
+ return origin;
+}
+
+//______________________________________________________
+Bool_t AliHFEdebugTreeTask::CheckITSstatus( const AliESDtrack * const esdtrack, Int_t layer) const {
+ //
+ // Check whether ITS area is dead
+ //
+ Int_t itsStatus = 0;
+ Int_t det;
+ Float_t xloc, zloc;
+ esdtrack->GetITSModuleIndexInfo(layer, det, itsStatus, xloc, zloc);
+ Bool_t status;
+ switch(itsStatus){
+ case 2: status = kFALSE; break;
+ case 3: status = kFALSE; break;
+ case 7: status = kFALSE; break;
+ default: status = kTRUE;
+ }
+ return status;
+}
+