#include <TObjArray.h>
#include <TClonesArray.h>
+#include <TPad.h>
#include <TProfile.h>
#include <TMath.h>
+#include <TDatabasePDG.h>
#include "TTreeStream.h"
#include "AliMagF.h"
#include "AliPID.h"
+#include "AliESDtrack.h"
#include "AliMathBase.h"
#include "AliTrackReference.h"
#include "AliAnalysisManager.h"
#include "AliTRDefficiencyMC.h"
ClassImp(AliTRDefficiencyMC)
-
+Float_t AliTRDefficiencyMC::fgPCut = 0.2; //[GeV/c]
+Float_t AliTRDefficiencyMC::fgPhiCut = 50.; //[deg]
+Float_t AliTRDefficiencyMC::fgThtCut = 50.; //[deg]
//_____________________________________________________________________________
AliTRDefficiencyMC::AliTRDefficiencyMC()
:AliTRDrecoTask("EfficiencyMC", "Combined Tracking Efficiency")
// Fill the histograms
//
const Int_t kArraySize = 10000; // Store indices of track references in an array
- Int_t indexAccepted[kArraySize], indexRejected[kArraySize], indexContamination[kArraySize];
- memset(indexAccepted, 0, sizeof(Int_t) * kArraySize);
- memset(indexRejected, 0, sizeof(Int_t) * kArraySize);
- memset(indexContamination, 0, sizeof(Int_t) * kArraySize);
- Int_t naccepted = 0, nrejected = 0, ncontamination = 0;
+ Int_t indexAccept[kArraySize],
+ indexReject[kArraySize],
+ indexContam[kArraySize];
+ memset(indexAccept, 0, sizeof(Int_t) * kArraySize);
+ memset(indexReject, 0, sizeof(Int_t) * kArraySize);
+ memset(indexContam, 0, sizeof(Int_t) * kArraySize);
+ Int_t naccept(0),
+ nreject(0),
+ nfindnt(0),
+ nkink(0),
+ ncontam(0);
Bool_t isContamination = kFALSE;
Int_t nTrackInfos = fTracks->GetEntriesFast();
- AliTRDtrackV1 *trackTRD = 0x0;
- AliTRDtrackInfo *trkInf = 0x0;
+ AliDebug(2, Form(" CANDIDATE TRACKS[%d]", nTrackInfos));
+
+ AliTRDtrackV1 *trackTRD(NULL);
+ AliTRDtrackInfo *trkInf(NULL);
for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
if(!trkInf) continue;
+
if(trkInf->GetTrack() || trkInf->GetNumberOfClustersRefit()){
- isContamination = IsRegistered(trkInf,indexAccepted,naccepted);
+ isContamination = (IsRegistered(trkInf,indexAccept,naccept)>=0);
if(!trkInf->GetNTrackRefs()){
// We reject the track since the Monte Carlo Information is missing
- AliDebug(1, Form("Error: Track Reference missing for Track %d", trkInf->GetLabel()));
+ AliDebug(2, Form("MC(Track Reference) missing @ label[%d]", trkInf->GetLabel()));
isContamination = kTRUE;
// Debugging
- if(trackTRD && DebugLevel() > 5) FillStreamTrackWOMC(trkInf);
- }
+ if(trackTRD && DebugLevel()>5) FillStreamTrackWOMC(trkInf);
+ }
if(isContamination){
// reject kink (we count these only once)
- if(trkInf->GetKinkIndex()) continue;
+ if(trkInf->GetKinkIndex()){
+ AliDebug(4, Form(" track @ idx[%d] MC[%d] is kink.", itinf, trkInf->GetLabel()));
+ nkink++;
+ continue;
+ }
// Register track as contamination
- indexContamination[ncontamination++]=itinf;
+ AliDebug(4, Form(" track @ idx[%d] MC[%d] is contamination.", itinf, trkInf->GetLabel()));
+ indexContam[ncontam++]=itinf;
continue;
}
// Accept track
- AliDebug(4, "Accept track");
+ AliDebug(4, Form(" track @ idx[%d] is ACCEPTED.", itinf));
// Register track as accepted
- indexAccepted[naccepted++] = itinf;
+ indexAccept[naccept++] = itinf;
}else{
- if(IsFindable(trkInf)){
+ Int_t code(0);
+ if((code=IsFindableNot(trkInf))){
+ AliDebug(4, Form(" track @ idx[%d] MC[%d] not findable [%d].", itinf, trkInf->GetLabel(), code));
+ nfindnt++;
+ } else {
// register track as rejected if not already registered there
// Attention:
// These track infos are not!!! registered as contamination
- if(!IsRegistered(trkInf, indexRejected, nrejected)) indexRejected[nrejected++] = itinf;
+ if(IsRegistered(trkInf, indexReject, nreject)<0){
+ AliDebug(4, Form(" track @ idx[%d] MC[%d] is missed.", itinf, trkInf->GetLabel()));
+ indexReject[nreject++] = itinf;
+ }
}
}
}
+ AliDebug(2, Form("TRACKS STATISTICS naccept[%d] ncontam[%d] nkink[%d] nmissed[%d] nfindnt[%d] ALL[%d] LOST[%d]", naccept, ncontam, nkink, nreject, nfindnt, naccept+ncontam+nkink, nTrackInfos-(naccept+nreject+ncontam+nkink+nfindnt)));
+
// we have to check if the rejected tracks are registered as found
// a possible source for this:
// first the track is not found by the barrel tracking but it is later found
// label would be created
// Attention:
// these tracks are not! registered as contamination
- Int_t tmprejected[kArraySize]; Int_t nrej = nrejected;
- memcpy(tmprejected, indexRejected, sizeof(Int_t) * nrejected);
- nrejected = 0;
+ Int_t tmprejected[kArraySize]; Int_t nrej = nreject;
+ memcpy(tmprejected, indexReject, sizeof(Int_t) * nreject);
+ nreject = 0;
for(Int_t irej = 0; irej < nrej; irej++){
trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej]));
- if(!IsRegistered(trkInf,indexAccepted,naccepted)) indexRejected[nrejected++] = tmprejected[irej];
+ Int_t idx(-1);
+ if((idx = IsRegistered(trkInf,indexAccept,naccept))<0){
+ indexReject[nreject++] = tmprejected[irej];
+ }else{
+ //printf("tracks @ accept[%d] missed[%d] are the same.\n", indexAccept[idx], tmprejected[irej]);
+ }
}
+
// Fill Histograms
- FillHistograms(naccepted, &indexAccepted[0], kAccepted);
- FillHistograms(nrejected, &indexRejected[0], kRejected);
- FillHistograms(ncontamination, &indexContamination[0], kContamination);
- Int_t nall = naccepted + nrejected;
- //if(DebugLevel()>=1)
- AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, nall ? 1.E2*Float_t(naccepted)/Float_t(nall) : 0.));
- AliInfo(Form("%3d Tracks: ALL[%3d] Contamination[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ncontamination, ncontamination, nall ? 1.E2*Float_t(ncontamination)/Float_t(nall + ncontamination) : 0.));
+ FillHistograms(naccept, &indexAccept[0], kAccept);
+ FillHistograms(nreject, &indexReject[0], kMiss);
+ FillHistograms(ncontam, &indexContam[0], kFake);
+
+ Int_t nall(naccept + nreject);
+ AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] Fake[%3d | %5.2f%%]",
+ (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(),
+ nall,
+ naccept,
+ (nall ? 1.E2*Float_t(naccept)/Float_t(nall) : 0.),
+ ncontam,
+ (nall ? 1.E2*Float_t(ncontam)/Float_t(nall) : 0.)));
PostData(0, fContainer);
}
//
// Change the histogram style
// For species histograms apply the colors connected with the given particle species
- //
- TH1 *histo = dynamic_cast<TH1 *>(fContainer->At(kEfficiencyHistogram));
- histo->SetMarkerStyle(22);
- histo->SetMarkerColor(kBlue);
- histo->GetXaxis()->SetTitle("p [GeV/c]");
- histo->GetXaxis()->SetMoreLogLabels();
- histo->GetYaxis()->SetTitle("Efficiency [%]");
- histo->GetYaxis()->SetRangeUser(0.99, 1.005);
-
- histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationHistogram));
- histo->SetMarkerStyle(22);
- histo->SetMarkerColor(kBlue);
- histo->GetXaxis()->SetTitle("p [GeV/c]");
- histo->GetXaxis()->SetMoreLogLabels();
- histo->GetYaxis()->SetTitle("Contamination [%]");
-
- // Species Efficiency Histograms
- for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
- histo = dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram + ispec));
- histo->SetMarkerStyle(22);
- histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
- histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
- histo->GetXaxis()->SetTitle("p [GeV/c]");
- histo->GetXaxis()->SetMoreLogLabels();
- histo->GetYaxis()->SetTitle("Efficiency [%]");
- histo->GetYaxis()->SetRangeUser(0.99, 1.005);
- }
-
- // Species Contamination Histograms
- for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
- histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram + ispec));
- histo->SetMarkerStyle(22);
- histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
- histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
- histo->GetXaxis()->SetTitle("p [GeV/c]");
- histo->GetXaxis()->SetMoreLogLabels();
- histo->GetYaxis()->SetTitle("Contamination [%]");
- }
-
- fNRefFigures = 6;
+ fNRefFigures = 8;
return kTRUE;
}
// Plot the histograms
//
if(ifig >= fNRefFigures) return kFALSE;
+ if(!gPad) return kFALSE;
+ gPad->SetLogx(kTRUE);
if(ifig < 2){
(dynamic_cast<TH1 *>(fContainer->At(ifig)))->Draw("e1");
return kTRUE;
switch(ifig){
case 2:
(dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram)))->Draw("e1");
- for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
- (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram + ispec)))->Draw("e1same");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+1)))->Draw("e1same");
break;
case 3:
- (dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram)))->Draw("e1");
- for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
- (dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram + ispec)))->Draw("e1same");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+2)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+3)))->Draw("e1same");
break;
case 4:
- (dynamic_cast<TH1 *>(fContainer->At(kEfficiencyNoPID)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+4)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+5)))->Draw("e1same");
break;
case 5:
- (dynamic_cast<TH1 *>(fContainer->At(kContaminationNoPID)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+6)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+7)))->Draw("e1same");
+ break;
+ case 6:
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+8)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+9)))->Draw("e1same");
+ break;
+ case 7:
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+10)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+11)))->Draw("e1same");
+ break;
+ case 8:
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+12)))->Draw("e1");
+ (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+13)))->Draw("e1same");
break;
}
return kTRUE;
//
// Create the histograms
//
- const Int_t nbins = 11;
if(fContainer) return fContainer;
- Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
-
- TString species[AliPID::kSPECIES] = {"Pions", "Muons", "Electrons", "Kaons", "Protons"};
- TString species_short[AliPID::kSPECIES] = {"Pi", "Mu", "El", "Ka", "Pr"};
-
- fContainer = new TObjArray();
- fContainer->AddAt(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyHistogram);
- fContainer->AddAt(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins), kContaminationHistogram);
- for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
- fContainer->AddAt(new TProfile(Form("trEffComb%s", species_short[ispec].Data()), Form("Combined Tracking Efficiency %s", species[ispec].Data()), nbins, xbins), kEfficiencySpeciesHistogram + ispec);
- for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
- fContainer->AddAt(new TProfile(Form("trContComb%s", species_short[ispec].Data()), Form("Combined Tracking Contamination %s", species[ispec].Data()), nbins, xbins), kContaminationSpeciesHistogram + ispec);
- fContainer->AddAt(new TProfile("trEffCombNoPID", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyNoPID);
- fContainer->AddAt(new TProfile("trContCombNoPID", "Combined Tracking Contamination", nbins, xbins), kContaminationNoPID);
+ const Int_t nbins = AliTRDCalPID::kNMom;
+ Float_t xbins[nbins+1] = {fgPCut, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
+ const Int_t marker[2][AliPID::kSPECIES+1] = {
+ {20, 21, 22, 23, 29, 2},
+ {24, 25, 26, 27, 30, 5}
+ };
+
+ fContainer = new TObjArray();fContainer->Expand(14);
+
+ TH1 *h(NULL);
+ fContainer->AddAt(h=new TProfile("hEff", "Tracking Efficiency ALL", nbins, xbins), kEfficiencyHistogram);
+ h->SetMarkerStyle(22);
+ h->SetMarkerColor(kBlue);
+ h->GetXaxis()->SetTitle("p [GeV/c]");
+ h->GetXaxis()->SetMoreLogLabels();
+ h->GetYaxis()->SetTitle("Efficiency");
+ h->GetYaxis()->SetRangeUser(0.2, 1.1);
+ fContainer->AddAt(h=new TProfile("hFake", "Fake Tracks", nbins, xbins), kContaminationHistogram);
+ h->SetMarkerStyle(22);
+ h->SetMarkerColor(kBlue);
+ h->GetXaxis()->SetTitle("p [GeV/c]");
+ h->GetXaxis()->SetMoreLogLabels();
+ h->GetYaxis()->SetTitle("Contamination");
+
+ Char_t sign[]={'+', '-'};
+ for(Int_t isign = 0; isign < 2; isign++){
+ for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
+ fContainer->AddAt(h=new TProfile(
+ Form("hEff_%s%c", AliPID::ParticleShortName(ispec), sign[isign]),
+ Form("Tracking Efficiency for %s%c", AliPID::ParticleName(ispec), sign[isign]), nbins, xbins),
+ kEfficiencySpeciesHistogram+ispec*2+isign);
+ h->SetMarkerStyle(marker[isign][ispec]);
+ h->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
+ h->SetMarkerColor(kBlack);
+ h->GetXaxis()->SetTitle("p [GeV/c]");
+ h->GetXaxis()->SetMoreLogLabels();
+ h->GetYaxis()->SetTitle("Efficiency");
+ h->GetYaxis()->SetRangeUser(0.2, 1.1);
+ }
+
+ fContainer->AddAt(h=new TProfile(Form("hEff_PID%c", sign[isign]), Form("Tracking Efficiency no PID %c", sign[isign]), nbins, xbins), kEfficiencySpeciesHistogram+AliPID::kSPECIES*2+isign);
+ h->SetMarkerStyle(marker[isign][AliPID::kSPECIES]);
+ h->SetMarkerColor(kBlack);h->SetLineColor(kBlack);
+ h->GetXaxis()->SetTitle("p [GeV/c]");
+ h->GetXaxis()->SetMoreLogLabels();
+ h->GetYaxis()->SetTitle("Efficiency");
+ h->GetYaxis()->SetRangeUser(0.2, 1.1);
+ }
return fContainer;
}
//_____________________________________________________________________________
-Bool_t AliTRDefficiencyMC::IsFindable(AliTRDtrackInfo * const trkInf){
+Int_t AliTRDefficiencyMC::IsFindableNot(AliTRDtrackInfo * const trkInf){
//
// Apply Cuts on the Monte Carlo track references
// return whether track is findable or not
//
- const Float_t kAlpha = 0.349065850;
- AliDebug(10, "Analysing Track References\n");
+
+ const Float_t chmbHght = AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
+ const Float_t eps(1.E-3);
+ Int_t ntr(trkInf->GetNTrackRefs());
+
+ AliDebug(10, Form(" CANDIDATE TrackRefs[%d]", ntr));
// Check if track is findable
- Double_t mom = 0.;
+ Double_t mom(0.), phi(0.), tht(0.);
Float_t xmin = 10000.0, xmax = 0.0;
Float_t ymin = 0.0, ymax = 0.0;
Float_t zmin = 0.0, zmax = 0.0;
Float_t lastx = 0.0, x = 0.0;
- Int_t nLayers = 0;
+ Int_t nLayers(0), ntrTRD(0);
Int_t sector[20];
- AliTrackReference *trackRef = 0x0;
- for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
- trackRef = trkInf->GetTrackRef(itr);
- AliDebug(10, Form("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z()));
+ AliTrackReference *trackRef(NULL);
+ for(Int_t itr(0); itr<ntr; itr++){
+ if(!(trackRef = trkInf->GetTrackRef(itr))) continue;
x = trackRef->LocalX();
-
// Be Sure that we are inside TRD
if(x < AliTRDinfoGen::GetTPCx() || x > AliTRDinfoGen::GetTOFx()) continue;
- sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
+ sector[ntrTRD] = Int_t(trackRef->Alpha()/AliTRDgeometry::GetAlpha());
+ AliDebug(10, Form(" [%2d] x[%7.2f] y[%7.2f] z[%7.2f] Sec[%2d]", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z(), sector[ntrTRD]));
if(x < xmin){
xmin = trackRef->LocalX();
ymin = trackRef->LocalY();
zmin = trackRef->Z();
- mom = trackRef->P();
+ mom = trackRef->P();
} else if(x > xmax){
xmax = trackRef->LocalX();
ymax = trackRef->LocalY();
}
if(itr > 0){
Float_t dist = TMath::Abs(x - lastx);
- AliDebug(10, Form("x = %f, lastx = %f, dist = %f\n", x, lastx, dist));
- if(TMath::Abs(dist - 3.7) < 0.1) nLayers++; // ref(i+1) has to be larger than ref(i)
+ if(TMath::Abs(dist - chmbHght) < eps && sector[ntrTRD]==sector[0]){
+ AliDebug(10, Form(" dx = %7.2f", dist));
+ nLayers++;
+ }
}
lastx = x;
+ ntrTRD++; if(ntrTRD>=20) break;
+ }
+ Double_t dx(xmax - xmin);
+ if(TMath::Abs(dx)<eps) return kNoChmb;
+
+ phi = (ymax -ymin)/dx;
+ tht = (zmax -zmin)/dx;
+ phi=TMath::ATan(phi)*TMath::RadToDeg();
+ tht=TMath::ATan(tht)*TMath::RadToDeg();
+ Bool_t primary = trkInf->IsPrimary();
+ const AliTRDtrackInfo::AliESDinfo *esd(trkInf->GetESDinfo());
+ AliDebug(10, Form(" p=%6.3f[GeV/c] phi=%6.2f[deg] theta=%6.2f[deg] nLy[%d]",
+ mom, phi, tht, nLayers));
+ if(DebugLevel()){
+ (*DebugStream()) << "IsFindable"
+ << "P=" << mom
+ << "Phi=" << phi
+ << "Tht=" << tht
+ << "Ntr=" << ntrTRD
+ << "NLy=" << nLayers
+ << "Primary=" << primary
+ << "\n";
}
// Apply cuts
- Bool_t findable = kTRUE;
- if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
- if(mom < 0.55) findable = kFALSE; // momentum cut at 0.6
- Double_t yangle = (ymax -ymin)/(xmax - xmin);
- Double_t zangle = (zmax -zmin)/(xmax - xmin);
- AliDebug(10, Form("\ntrack: y-Angle = %f, z-Angle = %f\nnLayers = %d", yangle, zangle, nLayers));
- if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
- if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
- if(nLayers < 4) findable = kFALSE;
- if(!trkInf->IsPrimary()) findable = kFALSE;
- Bool_t samesec = kTRUE;
- for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
- if(sector[iref] != sector[0]) samesec = kFALSE;
- if(!samesec) findable = kFALSE; // Discard all tracks which are passing more than one sector
- if(DebugLevel()){
- Double_t trackAngle = TMath::ATan(yangle);
- Bool_t primary = trkInf->IsPrimary();
- (*DebugStream()) << "NotFoundTrack"
- << "Momentum=" << mom
- << "trackAngle="<< trackAngle
- << "NLayers=" << nLayers
- << "Primary=" << primary
- << "\n";
- }
+ if(!nLayers) return kNoChmb;
+ if(xmax < xmin) return kCurved;
+ if(mom < fgPCut) return kPCut;
+
+
+ if(TMath::Abs(phi) > fgPhiCut) return kPhiCut;
+ if(TMath::Abs(tht) > fgThtCut) return kThtCut;
+
+ if(nLayers < 4){
+ if(!esd)return kLayer;
+ if(!(esd->GetStatus() & AliESDtrack::kTPCout)) return kLayer;
}
- else
- findable = kFALSE;
- return findable;
+
+ //if(!trkInf->IsPrimary()) {failCode=kPrimary; return kFALSE;}
+
+ return kFindable;
}
//_____________________________________________________________________________
-void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, FillingMode_t mode){
+void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, ETRDefficiencyMCstatus mode){
//
// Fill Histograms in three different modes:
// 1st tracks which are found and accepted
// 2nd tracks which are not found and not already counted
// 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
//
- const Int_t pid[AliPID::kSPECIES] = {211,13,11,321,2212};
- Double_t trkmom = 0.; // the track momentum
- Int_t trkpid = -1; // particle species
- AliTRDtrackInfo *trkInf = 0x0;
+
+ TDatabasePDG *dbPDG(TDatabasePDG::Instance());
+ Double_t trkmom(0.); // the track momentum
+ Int_t trkpdg(-1); // particle PDG code
+ AliTRDtrackInfo *trkInf(NULL);
for(Int_t itk = 0; itk < nTracks; itk++){
trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk]));
- AliDebug(3, Form("Accepted MC track: %d\n", trkInf->GetLabel()));
if(trkInf->GetNTrackRefs()){
// use Monte-Carlo Information for Momentum and PID
trkmom = trkInf->GetTrackRef(0)->P();
- trkpid = trkInf->GetPDG();
+ trkpdg = trkInf->GetPDG();
}else{
// Use TPC Momentum
trkmom = trkInf->GetTrack()->P();
}
+
+ const Char_t *cmode(NULL);
switch(mode){
- case kAccepted:
+ case kAccept:
(dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 1);
(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
+ cmode="ACCEPT";
break;
- case kRejected:
+ case kMiss:
(dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 0);
(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
+ cmode="MISS";
break;
- case kContamination:
+ case kFake:
(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 1);
+ cmode="FAKE";
break;
}
+ AliDebug(3, Form(" track[%d] MC[%d] Mode[%s]", indices[itk], trkInf->GetLabel(), cmode));
+
// Fill species histogram
- Int_t partSpec = -1;
- for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
- if(trkpid == pid[ispec]) partSpec = ispec;
- }
- if(partSpec >= 0){
- switch(mode){
- case kAccepted:
- (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 1);
- (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
- break;
- case kRejected:
- (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 0); (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
- break;
- case kContamination:
- (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 1);
- break;
- }
- } else {
- // The particle Type is not registered
- (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyNoPID)))->Fill(trkmom, 1);
- (dynamic_cast<TProfile *>(fContainer->At(kContaminationNoPID)))->Fill(trkmom, 1);
- }
+ Int_t idxSpec = AliTRDpidUtil::Pdg2Pid(TMath::Abs(trkpdg));
+ Int_t sign = dbPDG->GetParticle(trkpdg)->Charge() > 0. ? 1 : 0;
+ //printf("[%d]%s pdg[%d] sign[%d]\n", idxSpec, AliPID::ParticleName(idxSpec), trkpdg, sign);
+ if(idxSpec < 0) idxSpec = AliPID::kSPECIES;
+ (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + idxSpec*2+sign)))->Fill(trkmom, mode==kAccept?1:0);
}
}
Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
Int_t labelsTRD[180]; // Container for the cluster labels
Int_t sortlabels[360]; // Cluster Labels sorted according their occurancy
- AliTRDseedV1 *tracklet = 0x0;
- AliTRDcluster *c = 0x0;
- Int_t nclusters = 0x0;
+ AliTRDseedV1 *tracklet(NULL);
+ AliTRDcluster *c(NULL);
+ Int_t nclusters(0);
AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
tracklet = trackTRD->GetTracklet(il);
if(!tracklet) continue;
tracklet->ResetClusterIter();
- c = 0x0;
+ c = NULL;
while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
}
// Determine Label and Frequency
realPhi = realtrack->GetTrackRef(0)->Phi();
}
}
- (*DebugStream()) << "TrackingEffMCfake"
+ (*DebugStream()) << "EffMCfake"
<< "Event=" << event
<< "Label=" << label
<< "labelTRD=" << labelTRD
}
//_____________________________________________________________________________
-Bool_t AliTRDefficiencyMC::IsRegistered(AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
+Int_t AliTRDefficiencyMC::IsRegistered(AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
//
// Checks if track is registered in a given mode
//
- Bool_t isRegistered = kFALSE;
- for(Int_t il = 0; il < nTracks; il++){
- if((dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))->GetLabel() == trkInf->GetLabel()){
- isRegistered = kTRUE;
- break;
- }
+
+ Int_t label(trkInf->GetLabel());
+ for(Int_t il(nTracks); il--;){
+ if((dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))->GetLabel() == label) return il;
}
- return isRegistered;
+ return -1;
}
//____________________________________________________________________
AliTRDinfoGen::AliTRDinfoGen():
AliTRDrecoTask("infoGen", "MC-REC TRD-track list generator")
- ,fESDev(0x0)
- ,fMCev(0x0)
- ,fESDfriend(0x0)
- ,fTrackInfo(0x0)
- ,fEventInfo(0x0)
- ,fV0container(0x0)
- ,fV0Info(0x0)
+ ,fESDev(NULL)
+ ,fMCev(NULL)
+ ,fESDfriend(NULL)
+ ,fTrackInfo(NULL)
+ ,fEventInfo(NULL)
+ ,fV0container(NULL)
+ ,fV0Info(NULL)
{
//
// Default constructor
AliTRDinfoGen::~AliTRDinfoGen()
{
// Destructor
- if(fTrackInfo) delete fTrackInfo; fTrackInfo = 0x0;
- if(fEventInfo) delete fEventInfo; fEventInfo = 0x0;
- if(fV0Info) delete fV0Info; fV0Info = 0x0;
+
+ printf("AliTRDinfoGen::~AliTRDinfoGen() START\n");
+ if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;
+ printf("fTrackInfo\n");
+ if(fEventInfo) delete fEventInfo; fEventInfo = NULL;
+ printf("fEventInfo\n");
+ if(fV0Info) delete fV0Info; fV0Info = NULL;
+ printf("fV0Info\n");
if(fContainer){
fContainer->Delete(); delete fContainer;
- fContainer = 0x0;
+ fContainer = NULL;
}
+ printf("fContainer\n");
if(fV0container){
fV0container->Delete(); delete fV0container;
- fV0container = 0x0;
+ fV0container = NULL;
}
+ printf("fV0container\n");
+ printf("AliTRDinfoGen::~AliTRDinfoGen() STOP\n");
}
//____________________________________________________________________
//
TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
if(!tree){
- printf("ERROR - ESD event not found");
+ AliError("ESD event not found");
} else {
tree->SetBranchStatus("Tracks", 1);
tree->SetBranchStatus("ESDfriend*",1);
fESDev->SetESDfriend(fESDfriend);
new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
- Bool_t *trackMap = 0x0;
- AliStack * mStack = 0x0;
+ Bool_t *trackMap(NULL);
+ AliStack * mStack(NULL);
Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
if(HasMCdata()){
mStack = fMCev->Stack();
memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
}
- Int_t nTRD = 0, nTPC = 0, nclsTrklt;
+ Int_t nTRDout(0), nTRDin(0), nTPC(0), nclsTrklt;
AliDebug(2, Form("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC));
- AliESDtrack *esdTrack = 0x0;
- AliESDfriendTrack *esdFriendTrack = 0x0;
- TObject *calObject = 0x0;
- AliTRDtrackV1 *track = 0x0;
- AliTRDseedV1 *tracklet = 0x0;
- AliTRDcluster *cl = 0x0;
+ AliESDtrack *esdTrack = NULL;
+ AliESDfriendTrack *esdFriendTrack = NULL;
+ TObject *calObject = NULL;
+ AliTRDtrackV1 *track = NULL;
+ AliTRDseedV1 *tracklet = NULL;
+ AliTRDcluster *cl = NULL;
for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
esdTrack = fESDev->GetTrack(itrk);
AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));
- if(esdTrack->GetNcls(1)) nTPC++;
- if(esdTrack->GetNcls(2)) nTRD++;
+ if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
+ if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
+ if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;
// look at external track param
const AliExternalTrackParam *op = esdTrack->GetOuterParam();
AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
continue;
}
- AliMCParticle *mcParticle = 0x0;
+ AliMCParticle *mcParticle = NULL;
if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
continue;
}
fPdg = mcParticle->Particle()->GetPdgCode();
Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
- Int_t iref = 0; AliTrackReference *ref = 0x0;
+ Int_t iref = 0; AliTrackReference *ref = NULL;
while(iref<nRefs){
ref = mcParticle->GetTrackReference(iref);
if(ref->LocalX() > fgkTPC) break;
while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
- nTRD++;
AliDebug(4, Form("TRD track OK"));
// Set the clusters to unused
for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
fTrackInfo->Delete("");
}
- AliDebug(3, Form("%3d Tracks: TPC[%d] TRD[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRD));
+ AliDebug(2, Form("%3d Tracks: TPCout[%d] TRDin[%d] TRDout[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRDin, nTRDout));
-// AliESDv0 *v0 = 0x0;
+// AliESDv0 *v0 = NULL;
// for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
// if(!(v0 = fESD->GetV0(iv0))) continue;
// fV0container->Add(new AliTRDv0Info(v0));
AliMCParticle *mcParticle = (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));
Int_t fPdg = mcParticle->Particle()->GetPdgCode();
Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
- Int_t iref = 0; AliTrackReference *ref = 0x0;
+ Int_t iref = 0; AliTrackReference *ref = NULL;
Int_t nRefsTRD = 0;
new(fTrackInfo) AliTRDtrackInfo();
fTrackInfo->SetPDG(fPdg);
while(iref<nRefs){
+ Bool_t kIN(kFALSE);
ref = mcParticle->GetTrackReference(iref);
- if(ref->LocalX() > 250. && ref->LocalX() < 370.){
- AliDebug(4, Form(" trackRef[%2d] @ %7.3f IN", iref, ref->LocalX()));
+ if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTOF){
fTrackInfo->AddTrackRef(ref);
- nRefsTRD++;
+ nRefsTRD++;kIN=kTRUE;
}
- else AliDebug(4, Form(" trackRef[%2d] @ %7.3f OUT", iref, ref->LocalX()));
+ AliDebug(4, Form(" trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
iref++;
}
if(!nRefsTRD){
<< "\n";
info.Delete("");
}
- AliDebug(3, Form("Registering rejected MC track with label %d", itk));
+ AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
fContainer->Add(new AliTRDtrackInfo(*fTrackInfo));
fTrackInfo->Delete("");
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff