#include "AliLog.h"
#include "AliESDtrack.h"
#include "AliTrackReference.h"
+#include "AliAnalysisManager.h"
#include "AliTRDReconstructor.h"
#include "AliTRDtrackV1.h"
#include "AliTRDseedV1.h"
#include "AliTRDpidRefMaker.h"
+#include "AliTRDinfoGen.h"
+#include "info/AliTRDeventInfo.h"
#include "info/AliTRDv0Info.h"
+#include "info/AliTRDpidInfo.h"
// Defines and implements basic functionality for building reference data for TRD PID.
//
ClassImp(AliTRDpidRefMaker)
-ClassImp(AliTRDpidRefMaker::AliTRDpidRefData)
-ClassImp(AliTRDpidRefMaker::AliTRDpidRefDataArray)
-
-//________________________________________________________________________
-AliTRDpidRefMaker::AliTRDpidRefDataArray::AliTRDpidRefDataArray() :
- fNtracklets(0)
- ,fData(NULL)
-{
- // Constructor of data array
- fData = new AliTRDpidRefData[AliTRDgeometry::kNlayer];
-}
-
-//________________________________________________________________________
-AliTRDpidRefMaker::AliTRDpidRefDataArray::~AliTRDpidRefDataArray()
-{
- // Destructor
- delete [] fData;
-}
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::AliTRDpidRefDataArray::PushBack(Int_t ly, Int_t p, Float_t *dedx)
-{
-// Add PID data to the end of the array
- fData[fNtracklets].fPLbin= (ly<<4) | (p&0xf);
- memcpy(fData[fNtracklets].fdEdx, dedx, 8*sizeof(Float_t));
- fNtracklets++;
-}
-
-//________________________________________________________________________
-void AliTRDpidRefMaker::AliTRDpidRefDataArray::Reset()
-{
-// Reset content
-
- if(!fNtracklets) return;
- while(fNtracklets--){
- fData[fNtracklets].fPLbin = 0xff;
- memset(fData[fNtracklets].fdEdx, 0, 8*sizeof(Float_t));
- }
- fNtracklets=0;
-}
//________________________________________________________________________
AliTRDpidRefMaker::AliTRDpidRefMaker()
:AliTRDrecoTask()
- ,fReconstructor(NULL)
,fV0s(NULL)
,fData(NULL)
,fInfo(NULL)
,fPIDdataArray(NULL)
- ,fRefPID(kMC)
- ,fRefP(kMC)
- ,fPIDbin(0xff)
+ ,fRefPID(kV0)
+ ,fRefP(kRec)
,fFreq(1.)
,fP(-1.)
,fPthreshold(0.)
//________________________________________________________________________
AliTRDpidRefMaker::AliTRDpidRefMaker(const char *name, const char *title)
:AliTRDrecoTask(name, title)
- ,fReconstructor(NULL)
,fV0s(NULL)
,fData(NULL)
,fInfo(NULL)
,fPIDdataArray(NULL)
- ,fRefPID(kMC)
- ,fRefP(kMC)
- ,fPIDbin(0xff)
+ ,fRefPID(kV0)
+ ,fRefP(kRec)
,fFreq(1.)
,fP(-1.)
,fPthreshold(0.5)
// Default constructor
//
- fReconstructor = new AliTRDReconstructor();
- fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
- memset(fdEdx, 0, 10*sizeof(Float_t));
+ memset(fdEdx, 0, AliTRDpidUtil::kNNslices*sizeof(Float_t));
memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
- DefineInput(2, TObjArray::Class()); // v0 list
- DefineInput(3, TObjArray::Class()); // pid info list
+ DefineInput(3, TObjArray::Class()); // v0 list
+ DefineInput(4, TObjArray::Class()); // pid info list
DefineOutput(2, TTree::Class());
}
AliTRDpidRefMaker::~AliTRDpidRefMaker()
{
if(fPIDdataArray) delete fPIDdataArray;
- if(fReconstructor) delete fReconstructor;
}
//________________________________________________________________________
}
-//________________________________________________________________________
-void AliTRDpidRefMaker::ConnectInputData(Option_t *opt)
-{
- AliTRDrecoTask::ConnectInputData(opt);
- fV0s = dynamic_cast<TObjArray*>(GetInputData(2));
- fInfo = dynamic_cast<TObjArray*>(GetInputData(3));
-}
-
//________________________________________________________________________
void AliTRDpidRefMaker::UserCreateOutputObjects()
{
// Create histograms
// Called once
- OpenFile(1, "RECREATE");
fContainer = new TObjArray();
fContainer->SetName(Form("Moni%s", GetName()));
for(Int_t is=AliPID::kSPECIES; is--;) ax->SetBinLabel(is+1, AliPID::ParticleShortName(is));
h2->GetYaxis()->SetTitle("P bins");
h2->GetYaxis()->SetNdivisions(511);
-
fContainer->AddAt(h2, 0);
+ PostData(1, fContainer);
- fData = new TTree(GetName(), Form("Reference data for %s", GetName()));
- fData->Branch("s", &fPIDbin, "s/b");
+ OpenFile(2);
+ fData = new TTree("RefPID", "Reference data for PID");
fData->Branch("data", &fPIDdataArray);
+ PostData(2, fData);
}
//________________________________________________________________________
{
// Main loop
// Called for each event
+ Int_t ev((Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry());
+ if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))){
+ AliDebug(3, Form("Missing tracks container in ev %d", ev));
+ return;
+ }
+ if(!(fEvent = dynamic_cast<AliTRDeventInfo*>(GetInputData(2)))){
+ AliDebug(3, Form("Missing Event Info container in ev %d", ev));
+ return;
+ }
+ if(!(fV0s = dynamic_cast<TObjArray*>(GetInputData(3)))){
+ AliDebug(3, Form("Missing v0 container in ev %d", ev));
+ return;
+ }
+ if(!(fInfo = dynamic_cast<TObjArray*>(GetInputData(4)))){
+ AliDebug(3, Form("Missing pid info container in ev %d", ev));
+ return;
+ }
- AliInfo(Form("Analyse N[%d] tracks", fTracks->GetEntriesFast()));
-
+ AliDebug(1, Form("Entries: Ev[%d] Tracks[%d] V0[%d] PID[%d]", ev, fTracks->GetEntriesFast(), fV0s->GetEntriesFast(), fInfo->GetEntriesFast()));
AliTRDtrackInfo *track = NULL;
- AliTRDtrackV1 *trackTRD = NULL;
+ //AliTRDtrackV1 *trackTRD = NULL;
AliTrackReference *ref = NULL;
const AliTRDtrackInfo::AliESDinfo *infoESD = NULL;
for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
if(!track->HasESDtrack()) continue;
- trackTRD = track->GetTrack();
+ //trackTRD = track->GetTrack();
infoESD = track->GetESDinfo();
Double32_t *infoPID = infoESD->GetSliceIter();
Int_t n = infoESD->GetNSlices() - AliTRDgeometry::kNlayer;
AliDebug(4, Form("n[%d] p[GeV/c]{%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f}", n, p[0], p[1], p[2], p[3], p[4], p[5]));
ULong_t status = track->GetStatus();
- if(!(status&AliESDtrack::kTPCout)) continue;
+ if(!(status&AliESDtrack::kTRDpid)) continue;
// fill the pid information
- SetRefPID(fRefPID, track, fPID);
-
+ SetRefPID(fRefPID, track, infoESD, fPID);
+ // get particle type
+ Int_t idx(TMath::Max(Long64_t(0), TMath::LocMax(AliPID::kSPECIES, fPID)));
+ if(fPID[idx]<1.e-5) continue;
+
// prepare PID data array
if(!fPIDdataArray){
- fPIDdataArray = new AliTRDpidRefDataArray();
+ fPIDdataArray = new AliTRDpidInfo();
} else fPIDdataArray->Reset();
+ fPIDdataArray->SetPID(idx);
// fill PID information
for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
// fill P & dE/dx information
- if(HasFriends()){ // from TRD track
- if(!trackTRD) continue;
- AliTRDseedV1 *trackletTRD(NULL);
- trackTRD -> SetReconstructor(fReconstructor);
- if(!(trackletTRD = trackTRD -> GetTracklet(ily))) continue;
- if(!CheckQuality(trackletTRD)) continue;
- if(!CookdEdx(trackletTRD)) continue;
-
- // fill momentum information
- fP = 0.;
- switch(fRefP){
- case kMC:
- if(!(ref = track->GetTrackRef(trackletTRD))) continue;
- fP = ref->P();
- break;
- case kRec:
- fP = trackletTRD->GetMomentum();
- break;
- default: continue;
- }
- } else { // from ESD track
- // fill momentum information
- switch(fRefP){
- case kMC:
- if(!(ref = track->GetTrackRef(ily))) continue;
- fP = ref->P();
- break;
- case kRec:
- fP = p[ily];
- break;
- default: continue;
- }
- Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
- for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
+ switch(fRefP){
+ case kMC:
+ if(!(ref = track->GetTrackRef(ily))) continue;
+ fP = ref->P();
+ break;
+ case kRec:
+ fP = p[ily];
+ break;
+ default:
+ continue;
}
-
+ Double32_t *it = &infoPID[ily*AliTRDCalPID::kNSlicesNN];
+ for(Int_t is=AliTRDCalPID::kNSlicesNN; is--; it++) fdEdx[is] = (*it);
+
// momentum threshold
if(fP < fPthreshold) continue;
Fill();
}
-
- PostData(1, fContainer);
- PostData(2, fData);
}
{
// Fill data tree
- if(!fPIDdataArray->fNtracklets) return;
- fPIDbin = TMath::LocMax(AliPID::kSPECIES, fPID); // get particle type
-// Fill data tree
- AliInfo(Form("fData[%p]", (void*)fData));
+ if(!fPIDdataArray->GetNtracklets()) return;
+ // Fill data tree
fData->Fill();
// fill monitor
- for(Int_t ily=fPIDdataArray->fNtracklets; ily--;){
- Int_t pBin = fPIDdataArray->fData[ily].fPLbin & 0xf;
- ((TH2*)fContainer->At(0))->Fill(fPIDbin, pBin);
+ for(Int_t itrklt=fPIDdataArray->GetNtracklets(); itrklt--;){
+ Int_t pBin = fPIDdataArray->GetData(itrklt)->Momentum();
+ ((TH2*)fContainer->At(0))->Fill(fPIDdataArray->GetPID(), pBin);
}
}
void AliTRDpidRefMaker::LinkPIDdata()
{
// Link data tree to data members
- fData->SetBranchAddress("s", &fPIDbin);
fData->SetBranchAddress("data", &fPIDdataArray);
}
//________________________________________________________________________
-void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid)
+void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, AliTRDtrackInfo *track, const AliTRDtrackInfo::AliESDinfo *infoESD, Float_t *pid)
{
// Fill the reference PID values "pid" from "source" object
// according to the option "select". Possible options are
// - kMC - MC truth [default]
// - kRec - outside detectors
-
- memset(fPID, 0, AliPID::kSPECIES*sizeof(Float_t));
+ if(!track){
+ AliError("No trackInfo found");
+ return;
+ }
+ memset(pid, 0, AliPID::kSPECIES*sizeof(Float_t));
switch(select){
case kV0:
{
- AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
- if(!track){
- AliError("No trackInfo found");
- return;
- }
-
- //Get V0 PID decisions from the AliTRDv0Info for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) :
- AliTRDv0Info v0;
- for(Int_t is=AliPID::kSPECIES; is--;) fPID[is] = v0.GetV0PID(is, track);
+ //Get V0 PID decisions for all particle species (implemented so far : electrons from conversions, pions from K0s and protons from Lambdas) :
+ if(!infoESD->HasV0()) return;
+ const Int_t *v0pid=infoESD->GetV0pid();
+ for(Int_t is=AliPID::kSPECIES; is--;){ pid[is] = (Float_t)v0pid[is];}
}
break;
case kMC:
AliError("Could not retrive reference PID from MC");
return;
}
- {
- AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
- switch(track->GetPDG()){
- case kElectron:
- case kPositron:
- fPID[AliPID::kElectron] = 1.;
- break;
- case kMuonPlus:
- case kMuonMinus:
- fPID[AliPID::kMuon] = 1.;
- break;
- case kPiPlus:
- case kPiMinus:
- fPID[AliPID::kPion] = 1.;
- break;
- case kKPlus:
- case kKMinus:
- fPID[AliPID::kKaon] = 1.;
- break;
- case kProton:
- case kProtonBar:
- fPID[AliPID::kProton] = 1.;
- break;
- }
+ switch(track->GetPDG()){
+ case kElectron:
+ case kPositron:
+ pid[AliPID::kElectron] = 1.;
+ break;
+ case kMuonPlus:
+ case kMuonMinus:
+ pid[AliPID::kMuon] = 1.;
+ break;
+ case kPiPlus:
+ case kPiMinus:
+ pid[AliPID::kPion] = 1.;
+ break;
+ case kKPlus:
+ case kKMinus:
+ pid[AliPID::kKaon] = 1.;
+ break;
+ case kProton:
+ case kProtonBar:
+ pid[AliPID::kProton] = 1.;
+ break;
}
break;
case kRec:
{
- AliTRDtrackInfo *track = static_cast<AliTRDtrackInfo*>(source);
AliTRDtrackV1 *trackTRD = track->GetTrack();
- trackTRD -> SetReconstructor(fReconstructor);
+ trackTRD -> SetReconstructor(AliTRDinfoGen::Reconstructor());
//fReconstructor -> SetOption("nn");
trackTRD -> CookPID();
for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
AliWarning("PID reference source not implemented");
return;
}
- AliDebug(4, Form("Ref PID [%] : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]"
- ,AliPID::ParticleShortName(0), 1.e2*fPID[0]
- ,AliPID::ParticleShortName(1), 1.e2*fPID[1]
- ,AliPID::ParticleShortName(2), 1.e2*fPID[2]
- ,AliPID::ParticleShortName(3), 1.e2*fPID[3]
- ,AliPID::ParticleShortName(4), 1.e2*fPID[4]
+ AliDebug(4, Form("Ref PID : %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f] %s[%5.2f]"
+ ,AliPID::ParticleShortName(0), 1.e2*pid[0]
+ ,AliPID::ParticleShortName(1), 1.e2*pid[1]
+ ,AliPID::ParticleShortName(2), 1.e2*pid[2]
+ ,AliPID::ParticleShortName(3), 1.e2*pid[3]
+ ,AliPID::ParticleShortName(4), 1.e2*pid[4]
));
}