kNtrackletsFindable = 6,
kNtracksEvent = 7,
kNtracksSector = 8,
- kTrackStatus = 9,
- kTrackletStatus = 10,
- kPH = 11,
- kChi2 = 12,
- kChargeCluster = 13,
- kChargeTracklet = 14,
- kNeventsTrigger = 15,
- kNeventsTriggerTracks=16,
- kTriggerPurity = 17
+ kPH = 9,
+ kChi2 = 10,
+ kChargeCluster = 11,
+ kChargeTracklet = 12,
+ kNeventsTrigger = 13,
+ kNeventsTriggerTracks=14,
+ kTriggerPurity = 15,
+ kTrackStatus = 16,
+ kTrackletStatus = 17
};
AliTRDcheckDET();
}
//____________________________________________________________________
-TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t*)
+TGraphErrors* AliTRDcheckESD::GetGraph(Int_t id, Option_t *opt)
{
- Bool_t kBUILD = 1, // build graph if none found
- kCLEAR = 1; // clear existing graph
+ Bool_t kBUILD = strstr(opt, "b"), // build graph if none found
+ kCLEAR = strstr(opt, "c"); // clear existing graph
const Char_t *name[] = {
"Geo", "Trk", "Pid", "Ref"
TObjArray *res = 0x0;
if(!(res = (TObjArray*)fHistos->At(kResults)) ||
(id < 0 || id >= ngr)){
- AliWarning("Graph missing.");
+ AliWarning("Graph array missing.");
return 0x0;
}
h1[0] = h2->ProjectionX("px0", kTPCout, kTPCout);
h1[1] = h2->ProjectionX("px1", kTRDin, kTRDin);
Process(h1, GetGraph(0));
+ delete h1[0];delete h1[1];
// tracking efficiency
h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
h1[1] = h2->ProjectionX("px1", kTRDout, kTRDout);
Process(h1, GetGraph(1));
+ delete h1[1];
// PID efficiency
//h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
h1[1] = h2->ProjectionX("px1", kTRDpid, kTRDpid);
Process(h1, GetGraph(2));
+ delete h1[1];
// Refit efficiency
//h1[0] = h2->ProjectionX("px0", kTRDin, kTRDin);
h1[1] = h2->ProjectionX("px1", kTRDref, kTRDref);
Process(h1, GetGraph(3));
+ delete h1[1];
}
//____________________________________________________________________
void ConnectInputData(Option_t *);
void CreateOutputObjects();
- TGraphErrors* GetGraph(Int_t id, Option_t *opt=0x0);
+ TGraphErrors* GetGraph(Int_t id, Option_t *opt="bc");
void Exec(Option_t *);
Bool_t HasMC() const { return TESTBIT(fStatus, kMC);}
gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
Int_t bin = FindBin(is, 2.);
- h1 = h2->ProjectionY("px", bin, bin);
+ h1 = h2->ProjectionY(Form("px%d", is), bin, bin);
if(!h1->GetEntries()) continue;
h1->Scale(1./h1->Integral());
h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
FIRST = kTRUE;
for(Int_t is=0; is<AliPID::kSPECIES; is++){
Int_t bin = FindBin(is, 2.);
- h1 = h2->ProjectionY("pyt", bin, bin);
+ h1 = h2->ProjectionY(Form("pyt%d", is), bin, bin);
if(!h1->GetEntries()) continue;
h1->SetMarkerStyle(24);
h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
FIRST = kTRUE;
for(Int_t is=0; is<AliPID::kSPECIES; is++){
Int_t bin = FindBin(is, 2.);
- h1 = h2->ProjectionY("pyx", bin, bin);
+ h1 = h2->ProjectionY(Form("pyx%d", is), bin, bin);
if(!h1->GetEntries()) continue;
h1->SetMarkerStyle(24);
h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
legNClus->SetHeader("Particle Species");
for(Int_t is=0; is<AliPID::kSPECIES; is++){
Int_t bin = FindBin(is, 2.);
- h1 = h2->ProjectionY("pyNClus", bin, bin);
+ h1 = h2->ProjectionY(Form("pyNClus%d", is), bin, bin);
if(!h1->GetEntries()) continue;
h1->Scale(100./h1->Integral());
//h1->SetMarkerStyle(24);
legNClus->SetHeader("Particle Species");
for(Int_t is=0; is<AliPID::kSPECIES; is++){
Int_t bin = FindBin(is, 2.);
- h1 = h2->ProjectionY("pyNTracklets", bin, bin);
+ h1 = h2->ProjectionY(Form("pyNTracklets%d", is), bin, bin);
if(!h1->GetEntries()) continue;
h1->Scale(100./h1->Integral());
//h1->SetMarkerStyle(24);
}
if(fDebugLevel>=2) printf("%3d Tracks: TPC[%d] TRD[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRD);
- AliESDv0 *v0 = 0x0;
- for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
- if(!(v0 = fESD->GetV0(iv0))) continue;
- fV0container->Add(new AliTRDv0Info(v0));
- }
+// AliESDv0 *v0 = 0x0;
+// for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
+// if(!(v0 = fESD->GetV0(iv0))) continue;
+// fV0container->Add(new AliTRDv0Info(v0));
+// }
// Insert also MC tracks which are passing TRD where the track is not reconstructed
if(HasMCdata()){
//________________________________________________________________________
void AliTRDpidRefMaker::SetRefPID(ETRDpidRefMakerSource select, void *source, Float_t *pid)
{
-// !!!! PREMILMINARY FUNCTION !!!!
-//
-// this is the place for the V0 procedure
-// as long as there is no one implemented,
-// just the probabilities
-// of the TRDtrack are used!
+// Fill the reference PID values "pid" from "source" object
+// according to the option "select". Possible options are
+// - kV0 - v0 based PID
+// - kMC - MC truth [default]
+// - kRec - outside detectors
switch(select){
case kV0:
{
- AliTRDv0Info *v0 = static_cast<AliTRDv0Info*>(source);
- v0->Print(); // do something
+ 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);
}
break;
case kMC:
-
+#include "TMath.h"
+#include "AliESDtrack.h"
+#include "AliESDEvent.h"
#include "AliESDv0.h"
#include "AliTRDv0Info.h"
+#include "AliTRDtrackInfo.h"
+#include "AliESDInputHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliTRDtrackInfo.h"
+#include "AliLog.h"
+
+//Gathers all information necessary for reference data selection about the
+//track and (in case) its corresponding V0.
+//Carries out the selection of electrons (from gamma conversions), pions
+//(from K0s decays) and protons (from Lambda and Anti-Lambda decays) by
+//cuts specific for the respective decay and particle species.
+//(M.Heide, 2009/10/06)
+
+// Authors:
+// Alex Bercuci <A.Bercuci@gsi.de>
+// Alex Wilk <wilka@uni-muenster.de>
+// Markus Heide <mheide@uni-muenster.de>
+//
+
ClassImp(AliTRDv0Info)
//_________________________________________________
AliTRDv0Info::AliTRDv0Info()
: TObject()
- ,fStatus(0)
+ ,fESD(0x0)
+ ,fHasV0(0)
+ ,fQuality(0)
+
+ ,fMomentum(0)
+ ,fDCA(10)
+ ,fPointingAngle(10)
+ ,fOpenAngle(10)
+ ,fPsiPair(99)
+ ,fMagField(0)
+ ,fRadius(0)
+
+ ,fTrackID(0)
+ ,fV0Momentum(0)
+
+
+ ,fTrackP(0x0)
+ ,fTrackN(0x0)
+ ,fTrack(0x0)
+ ,fNindex(0)
+ ,fPindex(0)
+
{
+ memset(fPplus, 0, 2*kNlayer*sizeof(Float_t));
+ memset(fPminus, 0, 2*kNlayer*sizeof(Float_t));
+ memset(fDetPID, 0, 2*kNDaughters*kNDetectors*AliPID::kSPECIES*sizeof(Float_t));
+ memset(fInvMass, 0, kNMomBins*kNDecays*sizeof(Double_t));
+
+ /////////////////////////////////////////////////////////////////////////////
+ //Set Cut values: First specify decay in brackets, then the actual cut value!
+ /////////////////////////////////////////////////////////////////////////////
+
+ //Upper limit for distance of closest approach of two daughter tracks :
+ fUpDCA[kGamma] = 0.25;
+ fUpDCA[kK0s] = 0.25;
+ fUpDCA[kLambda] = 0.25;
+ fUpDCA[kAntiLambda] = 0.25;
+
+ //Upper limit for pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle) :
+ fUpPointingAngle[kGamma] = 0.03;
+ fUpPointingAngle[kK0s] = 0.03;
+ fUpPointingAngle[kLambda] = 0.03;
+ fUpPointingAngle[kAntiLambda] = 0.03;
+
+ //Upper limit for invariant mass of V0 mother :
+ fUpInvMass[kGamma][0] = 0.04;// second pair of brackets is for momentum bin: 0: below mother momentm of 2.5 GeV
+ fUpInvMass[kGamma][1] = 0.07;//1: above 2.5 GeV
+ fUpInvMass[kK0s][0] = fUpInvMass[kK0s][1] = 0.51;
+ fUpInvMass[kLambda][0] = fUpInvMass[kLambda][1] = 1.22;
+ fUpInvMass[kAntiLambda][0] = fUpInvMass[kAntiLambda][1] = 1.22;
+
+ //Lower limit for invariant mass of V0 mother :
+ fDownInvMass[kGamma] = -1.;
+ fDownInvMass[kK0s] = 0.49;
+ fDownInvMass[kLambda] = 1.;
+ fDownInvMass[kAntiLambda] = 1.;
+
+ //Lower limit for distance from secondary vertex to primary vertex in x-y plane :
+ fDownRadius[kGamma] = 5.7;
+ fDownRadius[kK0s] = 0.;
+ fDownRadius[kLambda] = 10.;
+ fDownRadius[kAntiLambda] = 10.;
+
+ //Upper limit for distance from secondary vertex to primary vertex in x-y plane :
+ fUpRadius[kGamma] = 1000.;
+ fUpRadius[kK0s] = 1000.;
+ fUpRadius[kLambda] = 1000.;
+ fUpRadius[kAntiLambda] = 1000.;
+
+ //Upper limit for opening angle between two daughter tracks (characteristically near zero for conversions) :
+ fUpOpenAngle[kGamma] = 0.1;
+ fUpOpenAngle[kK0s] = 3.15;
+ fUpOpenAngle[kLambda] = 3.15;
+ fUpOpenAngle[kAntiLambda] = 3.15;
+
+ //Upper limit for angle between daughter momentum plane and plane perpendicular to magnetic field (characteristically around zero for conversions) :
+ fUpPsiPair[kGamma] = 0.1;
+ fUpPsiPair[kK0s] = 1.6;
+ fUpPsiPair[kLambda] = 1.6;
+ fUpPsiPair[kAntiLambda] = 1.6;
+
+ //Lower limit for likelihood value of TPC PID :
+ fDownTPCPIDneg[AliPID::kElectron] = 0.21;
+ fDownTPCPIDpos[AliPID::kElectron] = 0.21;
+
+ fDownTPCPIDneg[AliPID::kMuon] = 0.21;
+ fDownTPCPIDpos[AliPID::kMuon] = 0.21;
+
+ fDownTPCPIDneg[AliPID::kPion] = 0.21;
+ fDownTPCPIDpos[AliPID::kPion] = 0.21;
+
+ fDownTPCPIDneg[AliPID::kKaon] = 0.21;
+ fDownTPCPIDpos[AliPID::kKaon] = 0.21;
+
+ fDownTPCPIDneg[AliPID::kProton] = 0.21;
+ fDownTPCPIDpos[AliPID::kProton] = 0.21;
+ //////////////////////////////////////////////////////////////////////////////////
+
}
//_________________________________________________
-AliTRDv0Info::AliTRDv0Info(AliESDv0 */*v0*/)
- : TObject()
- ,fStatus(0)
-{
+void AliTRDv0Info::GetESDv0Info(AliESDv0 *esdv0)
+{//Gets values of ESDv0 and daughter track properties
+ //See header file for description of variables
+
+ Int_t part1 = -1;
+ Int_t part2 = -1;
+
+ fQuality = Quality(esdv0);//Attributes an Int_t to the V0 due to quality cuts (= 1 if V0 is accepted, other integers depending on cut which excludes the vertex)
+
+ fRadius = Radius(esdv0);//distance from secondary vertex to primary vertex in x-y plane
+
+ fDCA = esdv0->GetDcaV0Daughters();//distance of closest approach of two daughter tracks
+
+ fPointingAngle = TMath::ACos(esdv0->GetV0CosineOfPointingAngle());// pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle)
+
+ fOpenAngle = OpenAngle(esdv0);//Opening angle between two daughter tracks
+
+ fPsiPair = PsiPair(esdv0);//Angle between daughter momentum plane and plane perpendicular to magnetic field
+
+ fV0Momentum = V0Momentum(esdv0);//Reconstructed momentum of the mother particle
+
+ for(Int_t idecay = 0; idecay < kNDecays; idecay++)//4 decay types : conversions, K0s, Lambda, Anti-Lambda
+ //five particle types: electrons, muons, pions, kaons, protons (muons and kaons not involved)
+ {
+ if(idecay == kLambda)//protons and pions from Lambda
+ {
+ part1 = AliPID::kProton;
+ part2 = AliPID::kPion;
+ }
+ else if(idecay == kAntiLambda)//antiprotons and pions from Anti-Lambda
+ {
+ part1 = AliPID::kPion;
+ part2 = AliPID::kProton;
+ }
+ else if(idecay == kK0s)//pions from K0s
+ part1 = part2 = AliPID::kPion;
+ else if(idecay == kGamma)//electrons from conversions
+ part1 = part2 = AliPID::kElectron;
+
+ fInvMass[idecay] = InvMass(part1, part2, esdv0);//Calculate invariant mass for all of our four supposed decays
+ }
+ GetDetectorPID();//Gets all likelihood values from TPC, TOF and ITS PID for the fDetPID[kNDaughters][kNDetectors][AliPID::kSPECIES] array
+
+
}
+//_________________________________________________
+Float_t AliTRDv0Info::V0Momentum(AliESDv0 *esdv0)
+{//Reconstructed momentum of V0 mother particle
+ Double_t mn[3] = {0,0,0};
+ Double_t mp[3] = {0,0,0};
+
+
+ esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+ esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
+
+
+ return TMath::Sqrt((mn[0]+mp[0])*(mn[0]+mp[0]) + (mn[1]+mp[1])*(mn[1]+mp[1])+(mn[2]+mp[2])*(mn[2]+mp[2]));
+}
+
+//_________________________________________________
+Double_t AliTRDv0Info::InvMass(Int_t part1, Int_t part2, AliESDv0 *esdv0)
+{//Invariant mass of reconstructed V0 mother
+
+ const Double_t kpmass[5] = {AliPID::ParticleMass(AliPID::kElectron),AliPID::ParticleMass(AliPID::kMuon),AliPID::ParticleMass(AliPID::kPion),AliPID::ParticleMass(AliPID::kKaon),AliPID::ParticleMass(AliPID::kProton)};
+ //Masses of electrons, muons, pions, kaons and protons, as implemented in ROOT
+
+
+ Double_t mn[3] = {0,0,0};
+ Double_t mp[3] = {0,0,0};
+
+ esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+ esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
+
+ Double_t mass1 = kpmass[part1];//sets supposed rest masses for both daughters
+ Double_t mass2 = kpmass[part2];
+
+ //Calculate daughters' energies :
+ Double_t e1 = TMath::Sqrt(mass1*mass1+
+ mp[0]*mp[0]+
+ mp[1]*mp[1]+
+ mp[2]*mp[2]);
+ Double_t e2 = TMath::Sqrt(mass2*mass2+
+ mn[0]*mn[0]+
+ mn[1]*mn[1]+
+ mn[2]*mn[2]);
+ //Sum of daughter momenta :
+ Double_t momsum =
+ (mn[0]+mp[0])*(mn[0]+mp[0])+
+ (mn[1]+mp[1])*(mn[1]+mp[1])+
+ (mn[2]+mp[2])*(mn[2]+mp[2]);
+ //invariant mass :
+ Double_t InvMass = TMath::Sqrt((e1+e2)*(e1+e2)-momsum);
+
+ return InvMass;
+
+}
+//_________________________________________________
+Float_t AliTRDv0Info::OpenAngle(AliESDv0 *esdv0)
+{//Opening angle between two daughter tracks
+ Double_t mn[3] = {0,0,0};
+ Double_t mp[3] = {0,0,0};
+
+
+ esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+ esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
+
+
+ fOpenAngle = TMath::ACos((mp[0]*mn[0] + mp[1]*mn[1] + mp[2]*mn[2])/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] + mp[2]*mp[2])*TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] + mn[2]*mn[2])));
+
+ return fOpenAngle;
+}
+
+//_________________________________________________
+Float_t AliTRDv0Info::PsiPair(AliESDv0 *esdv0)
+{//Angle between daughter momentum plane and plane perpendicular to magnetic field
+ Double_t x, y, z;
+ esdv0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
+
+ Double_t mn[3] = {0,0,0};
+ Double_t mp[3] = {0,0,0};
+
+
+ esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
+ esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
+
+
+ Double_t deltat = 1.;
+ deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
+
+ Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
+
+ Double_t MomPosProp[3];
+ Double_t MomNegProp[3];
+
+ AliExternalTrackParam nt(*fTrackN), pt(*fTrackP);
+
+ fPsiPair = 4.;
+
+ if(nt.PropagateTo(radiussum,fMagField) == 0)//propagate tracks to the outside
+ fPsiPair = -5.;
+ if(pt.PropagateTo(radiussum,fMagField) == 0)
+ fPsiPair = -5.;
+ pt.GetPxPyPz(MomPosProp);//Get momentum vectors of tracks after propagation
+ nt.GetPxPyPz(MomNegProp);
+
+ Double_t p_ele =
+ TMath::Sqrt(MomNegProp[0]*MomNegProp[0]+MomNegProp[1]*MomNegProp[1]+MomNegProp[2]*MomNegProp[2]);//absolute momentum value of negative daughter
+ Double_t p_pos =
+ TMath::Sqrt(MomPosProp[0]*MomPosProp[0]+MomPosProp[1]*MomPosProp[1]+MomPosProp[2]*MomPosProp[2]);//absolute momentum value of positive daughter
+
+ Double_t scalarproduct =
+ MomPosProp[0]*MomNegProp[0]+MomPosProp[1]*MomNegProp[1]+MomPosProp[2]*MomNegProp[2];//scalar product of propagated positive and negative daughters' momenta
+
+ Double_t chipair = TMath::ACos(scalarproduct/(p_ele*p_pos));//Angle between propagated daughter tracks
+
+ fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
+
+ return fPsiPair;
+
+}
+//_________________________________________________
+void AliTRDv0Info::V0fromTrack(AliTRDtrackInfo *track, Int_t ivertex)
+{//Checks if track is a secondary vertex daughter (due to V0 finder)
+
+ fMagField = fESD->GetMagneticField();
+
+ fTrackID = track->GetTrackId();//index of the track
+
+ fTrack = fESD->GetTrack(fTrackID);//sets track information
+
+ fHasV0 = 0;
+
+ //comparing index of track with indices of pos./neg. V0 daughter :
+ AliESDv0 * esdv0 = fESD->GetV0(ivertex);
+ if((esdv0->GetIndex(0) == fTrackID)||(esdv0->GetIndex(1) == fTrackID))
+ {
+ fHasV0 = 1;//track belongs to vertex found by V0 finder!
+ fNindex = esdv0->GetIndex(0);
+ fPindex = esdv0->GetIndex(1);
+ fTrackN = fESD->GetTrack(esdv0->GetIndex(0));//providing information about the other of the two daughter tracks
+ fTrackP = fESD->GetTrack(esdv0->GetIndex(1));
+ GetESDv0Info(esdv0);//gets all the relevant information about our V0
+ }
+}
+//_________________________________________________
+void AliTRDv0Info::GetDetectorPID()
+{//PID likelihoods from TPC, TOF, and ITS, for all particle species
+
+ fTrackN->GetTPCpid(fDetPID[kNeg][kTPC]);
+ fTrackP->GetTPCpid(fDetPID[kPos][kTPC]);
+ fTrackN->GetTOFpid(fDetPID[kNeg][kTOF]);
+ fTrackP->GetTOFpid(fDetPID[kPos][kTOF]);
+ fTrackN->GetITSpid(fDetPID[kNeg][kITS]);
+ fTrackP->GetITSpid(fDetPID[kPos][kITS]);
+
+}
+
+//_________________________________________________
+Float_t AliTRDv0Info::Radius(AliESDv0 *esdv0)
+{//distance from secondary vertex to primary vertex in x-y plane
+ Double_t x, y, z;
+ esdv0->GetXYZ(x,y,z); //Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
+ fRadius = TMath::Sqrt(x*x + y*y);
+ return fRadius;
+
+}
+//_________________________________________________
+Int_t AliTRDv0Info::Quality(AliESDv0 *esdv0)
+{ //Checking track and V0 quality status in order to exclude vertices based on poor information
+ Float_t NclsN;
+ NclsN = fTrackN->GetTPCNcls();//number of found clusters in TPC for negative track
+ Float_t NclsFN;
+ NclsFN = fTrackN->GetTPCNclsF();//number of findable clusters in TPC for negative track
+ Float_t NclsP;
+ NclsP = fTrackP->GetTPCNcls();//number of found clusters in TPC for positive track
+ Float_t NclsFP;
+ NclsFP = fTrackP->GetTPCNclsF();//number of findable clusters in TPC for positive track
+
+ fQuality = 0;
+
+
+ Float_t ClsRatioN;
+ Float_t ClsRatioP;
+
+ if((NclsFN == 0) || (NclsFP == 0))
+ return 2;
+
+ ClsRatioN = NclsN/NclsFN; //ratios of found to findable clusters in TPC
+ ClsRatioP = NclsP/NclsFP;
+
+ if (!(esdv0->GetOnFlyStatus()))//accept only vertices from online V0 finder
+ return 3;
+ if (!((fTrackP->GetStatus() &
+ AliESDtrack::kTPCrefit)))//accept only vertices in which both tracks have TPC refit
+ return 4;
+ if (!((fTrackN->GetStatus() &
+ AliESDtrack::kTPCrefit)))
+ return 5;
+ if (fTrackP->GetKinkIndex(0)>0 ||
+ fTrackN->GetKinkIndex(0)>0 )//exclude tracks with kinks
+ return 6;
+ if((ClsRatioN < 0.6)||(ClsRatioP < 0.6))//exclude tracks with low ratio of found to findable TPC clusters
+ return 7;
+ fQuality = 1;
+ return fQuality;
+}
+//_________________________________________________
+Bool_t AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track)
+{//decides if track is accepted for one of the reference data samples
+
+ Int_t iDecay = -1;
+
+ //decide which decay has to be considered for which particle sample (Anti-Lambda will be treated separately)
+ if(ipart == AliPID::kElectron)
+ iDecay = kGamma;
+ else if(ipart == AliPID::kPion)
+ iDecay = kK0s;
+ else if(ipart == AliPID::kProton)
+ iDecay = kLambda;
+
+ Int_t iPSlot;//Mother momentum slots above/below 2.5 GeV
+
+
+ Bool_t pid = 0;//return value for V0 PID decision
+
+ if(!(track))
+ {
+ AliError("AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track) : No track info found.\n");
+ return 0;
+ }
+
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ if(!esdH)
+ {
+ AliError("AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track) : ERROR - ESD input handler not found");
+ return 0;
+ }
+
+
+ fESD = esdH->GetEvent();
+
+ for(Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++)
+ {
+
+ if(pid == 0)
+ {
+ V0fromTrack(track, ivertex);//Get the V0 corresponding to the track (if there is a V0)
+
+ if(fV0Momentum > 2.5)//divide into slots according to reconstructed momentum of the mother particle
+ {iPSlot = 1;}
+ else
+ {iPSlot = 0;}
+ //Accept track for a sample only if...
+
+ if(!(fHasV0))//... there is a V0 found for it
+ continue;
+ if(!(fQuality == 1))//... it fulfills our quality criteria
+ continue;
+ if(!(fDCA < fUpDCA[iDecay]))//... distance of closest approach between daughters is reasonably small
+ continue;
+ if(!(fPointingAngle < fUpPointingAngle[iDecay]))//... pointing angle between momentum of mother particle and vector from prim. to sec. vertex is small
+ continue;
+ if(!(fRadius > fDownRadius[iDecay]))//... x-y plane distance of decay point to prim. vertex is bigger than a certain minimum value (for conversions)
+ continue;
+ if(!(fOpenAngle < fUpOpenAngle[iDecay]))//... opening angle is close enough to zero (for conversions)
+ continue;
+ if(!(TMath::Abs(fPsiPair) < fUpPsiPair[iDecay]))//... Psi-pair angle is close enough to zero(for conversions)
+ continue;
+
+ //specific cut criteria :
+ if(ipart == AliPID::kProton)
+ {//for proton sample: separate treatment of Lamba and Anti-Lambda decays:
+ //for Anti-Lambda:
+ //TPC PID likelihoods high enough for pi+ and anti-proton ; invariant mass calculated postulating these two particle species...
+ if((fDetPID[kNeg][kTPC][AliPID::kProton] > fDownTPCPIDneg[AliPID::kProton]) && (fDetPID[kPos][kTPC][AliPID::kPion] > fDownTPCPIDpos[AliPID::kPion]))
+ {
+ if(fNindex == fTrackID)
+ {
+ if((fInvMass[kAntiLambda] < fUpInvMass[kAntiLambda][iPSlot]) && (fInvMass[kAntiLambda] > fDownInvMass[kAntiLambda]))
+ {
+ pid = 1;
+ }
+ }
+ }
+ //for Lambda:
+ //TPC PID likelihoods high enough for pi- and proton ; invariant mass calculated accordingly
+ if((fDetPID[kNeg][kTPC][AliPID::kPion] > fDownTPCPIDneg[AliPID::kPion]) && (fDetPID[kPos][kTPC][AliPID::kProton] > fDownTPCPIDpos[AliPID::kProton]))
+ {
+ if(fPindex == fTrackID)
+ {
+ if((fInvMass[kLambda] < fUpInvMass[kLambda][iPSlot]) && (fInvMass[kLambda] > fDownInvMass[kLambda]))
+ {
+ pid = 1;
+ }
+ }
+ }
+ }
+ //for photon and K0s decays: equal TPC PID likelihood criteria for both daughters ; invariant mass calculated postulating two electrons/two pions
+ else
+ if((fDetPID[kNeg][kTPC][ipart] > fDownTPCPIDneg[ipart]) && (fDetPID[kPos][kTPC][ipart] > fDownTPCPIDpos[ipart]))
+ {
+ if((fInvMass[iDecay] < fUpInvMass[iDecay][iPSlot]) && (fInvMass[iDecay] > fDownInvMass[iDecay]))
+ {
+ pid = 1;
+ }
+ }
+ }
+ }
+
+ return pid;
+
+}
//_________________________________________________
void AliTRDv0Info::Print(Option_t */*opt*/) const
{
- printf("AliTRDv0Info::Print() : Nothing implemented so far.\n");
+
}
#include "AliPID.h"
#endif
+
class AliESDv0;
+class AliESDtrack;
+class AliESDEvent;
+class AliTRDtrackV1;
+class AliTRDtrackInfo;
class AliTRDv0Info : public TObject
{
public:
enum ETRDv0Info{
kNV0param = 10
- ,kNlayer = AliTRDgeometry::kNlayer
- ,kNDetectors = 2
+ ,kNlayer = AliTRDgeometry::kNlayer
+ ,kNDetectors = 3//TPC, TOF, ITS (TOF and ITS not implemented yet)
+ ,kNDaughters = 2//for positive and negative track
+ ,kNDecays = 4//number of decay types considered for reference data (conversions, K0s, Lambda, Anti-Lambda)
+ ,kNMomBins = 2//number of different momentum bins to consider for different cuts; first example: below/above 2.5 GeV -> to be refined!
+ };
+
+ enum EDaughter{
+ kNeg = 0
+ ,kPos = 1
+ };
+
+ enum EDecayType{
+ kGamma = 0
+ ,kK0s = 1
+ ,kLambda = 2
+ ,kAntiLambda = 3
};
+
+ enum EDetector{
+ kTPC = 0
+ ,kTOF = 1
+ ,kITS = 2
+ };
+
+
AliTRDv0Info();
- AliTRDv0Info(AliESDv0 *v0);
- void Print(Option_t *opt=0x0) const;
+ virtual ~AliTRDv0Info(){}
+ Float_t Pplus[2*kNlayer];
+ Float_t Pminus[2*kNlayer];
+
+
+ void Print(Option_t *opt=0x0) const;
+
+ Bool_t GetV0PID(Int_t ipart, AliTRDtrackInfo *track);//decides if a track is accepted for one of the reference samples!!
+
+ //Set values of measured/calculated variables:
+ void SetQuality(Int_t Quality){fQuality = Quality;}
+ void SetPplus(Int_t iLayer, Float_t Pplus){fPplus[iLayer] = Pplus;}
+ void SetPminus(Int_t iLayer, Float_t Pminus){fPminus[iLayer] = Pminus;}
+ void SetDCA(Float_t DCA){fDCA = DCA;}
+ void SetMomentum(Float_t Momentum){fMomentum = Momentum;}
+ void SetPointingAngle(Float_t PointingAngle){fPointingAngle = PointingAngle;}
+ void SetOpenAngle(Float_t OpenAngle){fOpenAngle = OpenAngle;}
+ void SetPsiPair(Float_t PsiPair){fPsiPair = PsiPair;}
+ void SetRadius(Float_t Radius){fRadius = Radius;}
+ void SetInvMass(Int_t iDecay, Float_t InvMass){fInvMass[iDecay] = InvMass;}
+ void SetDetPID(Int_t iDaughter, Int_t iDetector, Int_t iSpecies, Float_t DetPID){fDetPID[iDaughter][iDetector][iSpecies] = DetPID;}
+
+//____________________________________________________________
+ //Set cut values:
+
+ void SetUpDCA(Int_t iDecay, Float_t UpDCA){fUpDCA[iDecay] = UpDCA;}
+ void SetUpPointingAngle(Int_t iDecay, Float_t UpPointingAngle){fUpPointingAngle[iDecay] = UpPointingAngle;}
+ void SetUpOpenAngle(Int_t iDecay, Float_t UpOpenAngle){fUpOpenAngle[iDecay] = UpOpenAngle;}
+ void SetDownOpenAngle(Int_t iDecay, Float_t DownOpenAngle){fDownOpenAngle[iDecay] = DownOpenAngle;}
+ void SetUpPsiPair(Int_t iDecay, Float_t UpPsiPair){fUpPsiPair[iDecay] = UpPsiPair;}
+ void SetDownPsiPair(Int_t iDecay, Float_t DownPsiPair){fDownPsiPair[iDecay] = DownPsiPair;}
+ void SetUpRadius(Int_t iDecay, Float_t UpRadius){fUpRadius[iDecay] = UpRadius;}
+ void SetDownRadius(Int_t iDecay, Float_t DownRadius){fDownRadius[iDecay] = DownRadius;}
+ void SetUpInvMass(Int_t iDecay, Int_t iMomentum, Double_t UpInvMass){fUpInvMass[iDecay][iMomentum] = UpInvMass;}
+ void SetDownInvMass(Int_t iDecay, Double_t DownInvMass){fDownInvMass[iDecay] = DownInvMass;}
+ void SetDownTPCPIDneg(Int_t iDecay, Double_t DownTPCPIDneg){fDownTPCPIDneg[iDecay] = DownTPCPIDneg;}
+ void SetDownTPCPIDpos(Int_t iDecay, Double_t DownTPCPIDpos){fDownTPCPIDpos[iDecay] = DownTPCPIDpos;}
+
+
private:
- Int_t fStatus; // track status
- Float_t fV0param[kNV0param]; // V0 parameters (momentum, variance, etc)
+ AliTRDv0Info(const AliTRDv0Info&);
+ AliTRDv0Info& operator=(const AliTRDv0Info&);
+
+ void GetESDv0Info(AliESDv0 *esdv0);//gets most of the variables below
+ void GetDetectorPID();//operating with likelihood values of different detectors
+ Int_t Quality(AliESDv0 *esdv0);//checks for track/vertex quality criteria
+ Double_t InvMass(Int_t part1, Int_t part2, AliESDv0 *esdv0);//invariant mass of mother
+ Float_t PsiPair(AliESDv0 *esdv0);//angle between daughters in plane perpendicular to magnetic field (characteristically around zero for conversions)
+ Float_t OpenAngle(AliESDv0 *esdv0);//opening angle between V0 daughters; close to zero for conversions
+ Float_t Radius(AliESDv0 *esdv0);//distance of secondary to primary vertex in x-y-plane
+ Float_t DCA() const {return fDCA;}//distance of closest approach between supposed daughter tracks
+ Float_t PointingAngle() const {return fPointingAngle;}//pointing angle: between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle
+ Float_t V0Momentum(AliESDv0 *esdv0);//reconstructed momentum of V0 mother particle
+ void V0fromTrack(AliTRDtrackInfo *track, Int_t ivertex);//checks if a track belongs to a vertex found by V0 finder
+
+ AliESDEvent *fESD;
+
+
+ Bool_t fHasV0; //Does this track belong to a vertex from a V0 finder?
+
+ Int_t fQuality; // track quality status for both V0 daughters; OnFly, TPCrefit, Kinks, TPC clusters
+
Float_t fPplus[2*kNlayer]; // momentum and variance for the positive daughter
Float_t fPminus[2*kNlayer]; // momentum and variance for the negative daughter
- Float_t fPID[kNDetectors][AliPID::kSPECIES]; // PID provided by TPC and TOF
+ Double_t fDetPID[kNDaughters][kNDetectors][AliPID::kSPECIES]; // PID provided by TPC, TOF and ITS
+
+ Float_t fMomentum; // Momentum of track at the vertex
+
+ Float_t fDCA; // Distance of closest approach of daughter tracks
+
+ Float_t fPointingAngle;// = TMath::ACos(esdv0->GetV0CosineOfPointingAngle()); // Cosine of pointing angle
+
+ Float_t fOpenAngle; // opening angle between daughters
+
+ Float_t fPsiPair; // /Angle between daughter momentum plane and plane perpendicular to magnetic field
+
+ Double_t fInvMass[kNDecays]; // invariant mass for different decay scenarios (conversions, K0s, Lambda->p+pi-, Lambda->p-pi+)
+
+ Double_t fMagField;
+
+ Float_t fRadius; //distance of decay/conversion from primary vertex in x-y plane
+
+ Int_t fTrackID;//track index
+
+
+ Float_t fV0Momentum; //V0 mother's momentum
+
+ //____________________________________________________________
+ //Upper and lower limits for cut variables:
+
+ Float_t fUpDCA[kNDecays];
+
+ Float_t fUpPointingAngle[kNDecays];
+
+ Float_t fUpOpenAngle[kNDecays];
+
+ Float_t fDownOpenAngle[kNDecays];
+
+ Float_t fUpPsiPair[kNDecays];
+
+ Float_t fDownPsiPair[kNDecays];
+
+ Double_t fUpInvMass[kNDecays][kNMomBins];
+
+ Double_t fDownInvMass[kNDecays];
+
+ Float_t fUpRadius[kNDecays];
+
+ Float_t fDownRadius[kNDecays];
+
+ Float_t fDownTPCPIDneg[AliPID::kSPECIES];
+
+ Float_t fDownTPCPIDpos[AliPID::kSPECIES];
+
+
+ AliESDtrack *fTrackP; //positive daughter
+ AliESDtrack *fTrackN; //negative daughter
+ AliESDtrack *fTrack; //the current track in the ESDtrack loop (either positive or negative)
+
+ Int_t fNindex;//indices of positive and negative daughter track
+ Int_t fPindex;
+
+
ClassDef(AliTRDv0Info, 0) // extracted V0 MC information
};
void AddTRDcheckPID(AliAnalysisManager *mgr, Char_t *trd, AliAnalysisDataContainer **ci/*, AliAnalysisDataContainer **co*/)
{
Int_t map = ParseOptions(trd);
- if(!(TSTBIT(map, kCheckPID))) return;
+ if(TSTBIT(map, kCheckPID)){
+ AliTRDcheckPID *pid = 0x0;
+ mgr->AddTask(pid = new AliTRDcheckPID());
+ pid->SetDebugLevel(0);
+ pid->SetMCdata(mgr->GetMCtruthEventHandler());
+ mgr->ConnectInput(pid, 0, ci[0]);
+ mgr->ConnectOutput(pid, 0, mgr->CreateContainer(pid->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, "TRD.Performance.root"));
+ }
- AliTRDcheckPID *pid = 0x0;
- mgr->AddTask(pid = new AliTRDcheckPID());
- pid->SetDebugLevel(0);
- pid->SetMCdata(mgr->GetMCtruthEventHandler());
- mgr->ConnectInput(pid, 0, ci[0]);
- mgr->ConnectOutput(pid, 0, mgr->CreateContainer(pid->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, "TRD.Performance.root"));
-
- if(!(TSTBIT(map, kPIDRefMaker))) return;
- // TRD pid reference
- AliTRDpidRefMakerNN *ref = new AliTRDpidRefMakerNN();
- //AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 4);
- // Neural network PID
- mgr->AddTask(ref);
- ref->SetDebugLevel(3);
- AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 3);
- ref->SetMCdata(mgr->GetMCtruthEventHandler());
- mgr->ConnectInput( ref, 0, ci[0]);
- mgr->ConnectInput( ref, 1, ci[2]);
- mgr->ConnectOutput(ref, 0, mgr->CreateContainer(Form("Moni%s", ref->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
- mgr->ConnectOutput(ref, 1, mgr->CreateContainer(ref->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
+ if(TSTBIT(map, kPIDRefMakerNN)){
+ // TRD NN pid reference maker
+ AliTRDpidRefMakerNN *ref = new AliTRDpidRefMakerNN();
+ //AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 4);
+ // Neural network PID
+ mgr->AddTask(ref);
+ ref->SetDebugLevel(3);
+ AliLog::SetClassDebugLevel("AliTRDpidRefMakerNN", 3);
+ ref->SetMCdata(mgr->GetMCtruthEventHandler());
+ mgr->ConnectInput( ref, 0, ci[0]);
+ mgr->ConnectInput( ref, 1, ci[2]);
+ mgr->ConnectOutput(ref, 0, mgr->CreateContainer(Form("Moni%s", ref->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
+ mgr->ConnectOutput(ref, 1, mgr->CreateContainer(ref->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", ref->GetName())));
+ }
// Multidimensional Likelihood PID
- AliTRDpidRefMakerLQ *reflq = new AliTRDpidRefMakerLQ();
- mgr->AddTask(reflq);
- reflq->SetDebugLevel(3);
- //AliLog::SetClassDebugLevel("AliTRDpidRefMakerLQ", 3);
- reflq->SetMCdata(mgr->GetMCtruthEventHandler());
- mgr->ConnectInput( reflq, 0, ci[0]);
- mgr->ConnectInput( reflq, 1, ci[2]);
- mgr->ConnectOutput(reflq, 0, mgr->CreateContainer(Form("Moni%s", reflq->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
- mgr->ConnectOutput(reflq, 1, mgr->CreateContainer(reflq->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
+ if(TSTBIT(map, kPIDRefMakerLQ)){
+ AliTRDpidRefMakerLQ *reflq = new AliTRDpidRefMakerLQ();
+ mgr->AddTask(reflq);
+ reflq->SetDebugLevel(3);
+ //AliLog::SetClassDebugLevel("AliTRDpidRefMakerLQ", 3);
+ reflq->SetMCdata(mgr->GetMCtruthEventHandler());
+ mgr->ConnectInput( reflq, 0, ci[0]);
+ mgr->ConnectInput( reflq, 1, ci[2]);
+ mgr->ConnectOutput(reflq, 0, mgr->CreateContainer(Form("Moni%s", reflq->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
+ mgr->ConnectOutput(reflq, 1, mgr->CreateContainer(reflq->GetName(), TTree::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Calib%s.root", reflq->GetName())));
+ }
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff