Added corrected definition of cos(theta*) and related computation.
Corrected the cut for perfect PID.
#include "AliLog.h"
#include "AliVEvent.h"
+#include "AliMCEvent.h"
#include "AliRsnEvent.h"
#include "AliRsnPair.h"
#include "AliRsnAnalysisManager.h"
AliDebug(AliLog::kDebug+2,"->");
}
+
+//_____________________________________________________________________________
+void AliRsnAnalysisManager::ProcessAllPairsMC(AliRsnEvent *ev0, AliRsnEvent *ev1)
+{
+//
+// Process one or two events for all pair managers.
+//
+
+ AliDebug(AliLog::kDebug+2,"<-");
+
+ if (!ev1) ev1 = ev0;
+
+ Int_t nTracks[2];
+ nTracks[0] = ev0->GetRefMC()->GetNumberOfTracks();
+ nTracks[1] = ev1->GetRefMC()->GetNumberOfTracks();
+
+ // external loop
+ // joins the loop on tracks and v0s, by looping the indexes from 0
+ // to the sum of them, and checking what to take depending of its value
+ Int_t i0, i1, i;
+ AliRsnDaughter daughter0, daughter1;
+ AliRsnPair *pair = 0x0;
+ TObjArrayIter next(&fPairs);
+
+ for (i0 = 0; i0 < nTracks[0]; i0++)
+ {
+ // assign first track
+ if (i0 < nTracks[0]) ev0->SetDaughter(daughter0, i0, AliRsnDaughter::kTrack);
+ else ev0->SetDaughter(daughter0, i0 - nTracks[0], AliRsnDaughter::kV0);
+
+ // internal loop (same criterion)
+ for (i1 = 0; i1 < nTracks[1]; i1++)
+ {
+ // if looking same event, skip the case when the two indexes are equal
+ if (ev0 == ev1 && i0 == i1) continue;
+
+ // assign second track
+ if (i1 < nTracks[1]) ev1->SetDaughter(daughter1, i1, AliRsnDaughter::kTrack);
+ else ev1->SetDaughter(daughter1, i1 - nTracks[1], AliRsnDaughter::kV0);
+
+ // loop over all pairs and make computations
+ next.Reset();
+ i = 0;
+ while ((pair = (AliRsnPair*)next()))
+ {
+ AliDebug(AliLog::kDebug+1, Form("ProcessAllPairs of the AnalysisManager(%s) [%d] ...", pair->GetName(), i++));
+
+ // if the pair is a like-sign, skip the case when i1 < i0,
+ // in order not to double count each like-sign pair
+ // (equivalent to looping from i0+1 to ntracks)
+ if (pair->GetPairDef()->IsLikeSign() && i1 < i0) continue;
+
+ // process the two tracks
+ if (!pair->Fill(&daughter0, &daughter1, ev0, ev1)) continue;
+ pair->Compute();
+ }
+ }
+ }
+
+ AliDebug(AliLog::kDebug+2,"->");
+}
+
void InitAllPairs(TList*list);
void ProcessAllPairs(AliRsnEvent *ev0, AliRsnEvent *ev1);
+ void ProcessAllPairsMC(AliRsnEvent *ev0, AliRsnEvent *ev1);
private:
// the virtual class has already sorted tracks in the PID index
// so we need here just to call the execution of analysis
- fRsnAnalysisManager.ProcessAllPairs(&fRsnEvent, &fRsnEvent);
+ if (!fMCOnly) fRsnAnalysisManager.ProcessAllPairs(&fRsnEvent, &fRsnEvent);
+ else fRsnAnalysisManager.ProcessAllPairsMC(&fRsnEvent, &fRsnEvent);
PostData(2, fOutList);
AliDebug(AliLog::kDebug+2,"->");
fCutResult = (fCutValueI == fMinI);
// print debug message
- AliDebug(AliLog::kDebug + 3, "=== CUT DEBUG ====================================");
- AliDebug(AliLog::kDebug + 3, Form("Cut name : %s", GetName()));
- AliDebug(AliLog::kDebug + 3, Form("Checked value: %d", fCutValueI));
- AliDebug(AliLog::kDebug + 3, Form("Cut value : %d", fMinI));
- AliDebug(AliLog::kDebug + 3, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
- AliDebug(AliLog::kDebug + 3, "=== END CUT DEBUG ================================");
+ AliDebug(AliLog::kDebug + 2, "=== CUT DEBUG ====================================");
+ AliDebug(AliLog::kDebug + 2, Form("Cut name : %s", GetName()));
+ AliDebug(AliLog::kDebug + 2, Form("Checked value: %d", fCutValueI));
+ AliDebug(AliLog::kDebug + 2, Form("Cut value : %d", fMinI));
+ AliDebug(AliLog::kDebug + 2, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
+ AliDebug(AliLog::kDebug + 2, "=== END CUT DEBUG ================================");
return fCutResult;
}
fCutResult = (TMath::Abs(fCutValueD - fMinD) < 1E-6);
// print debug message
- AliDebug(AliLog::kDebug + 3, "=== CUT DEBUG ====================================");
- AliDebug(AliLog::kDebug + 3, Form("Cut name : %s", GetName()));
- AliDebug(AliLog::kDebug + 3, Form("Checked value: %f", fCutValueD));
- AliDebug(AliLog::kDebug + 3, Form("Cut value : %f", fMinD));
- AliDebug(AliLog::kDebug + 3, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
- AliDebug(AliLog::kDebug + 3, "=== END CUT DEBUG ================================");
+ AliDebug(AliLog::kDebug + 2, "=== CUT DEBUG ====================================");
+ AliDebug(AliLog::kDebug + 2, Form("Cut name : %s", GetName()));
+ AliDebug(AliLog::kDebug + 2, Form("Checked value: %f", fCutValueD));
+ AliDebug(AliLog::kDebug + 2, Form("Cut value : %f", fMinD));
+ AliDebug(AliLog::kDebug + 2, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
+ AliDebug(AliLog::kDebug + 2, "=== END CUT DEBUG ================================");
return fCutResult;
}
fCutResult = ((fCutValueI >= fMinI) && (fCutValueI <= fMaxI));
// print debug message
- AliDebug(AliLog::kDebug + 3, "=== CUT DEBUG ====================================");
+ AliDebug(AliLog::kDebug + 2, "=== CUT DEBUG ====================================");
AliDebug(AliLog::kDebug + 2, Form("Cut name : %s", GetName()));
AliDebug(AliLog::kDebug + 2, Form("Checked value: %d", fCutValueI));
AliDebug(AliLog::kDebug + 2, Form("Cut range : %d , %d", fMinI, fMaxI));
AliDebug(AliLog::kDebug + 2, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
- AliDebug(AliLog::kDebug + 3, "=== END CUT DEBUG ================================");
+ AliDebug(AliLog::kDebug + 2, "=== END CUT DEBUG ================================");
return fCutResult;
}
fCutResult = ((fCutValueD >= fMinD) && (fCutValueD <= fMaxD));
// print debug message
- AliDebug(AliLog::kDebug + 3, "=== CUT DEBUG ====================================");
+ AliDebug(AliLog::kDebug + 2, "=== CUT DEBUG ====================================");
AliDebug(AliLog::kDebug + 2, Form("Cut name : %s", GetName()));
AliDebug(AliLog::kDebug + 2, Form("Checked value: %f", fCutValueD));
AliDebug(AliLog::kDebug + 2, Form("Cut range : %f , %f", fMinD, fMaxD));
AliDebug(AliLog::kDebug + 2, Form("Cut result : %s", (fCutResult ? "PASSED" : "NOT PASSED")));
- AliDebug(AliLog::kDebug + 3, "=== END CUT DEBUG ================================");
+ AliDebug(AliLog::kDebug + 2, "=== END CUT DEBUG ================================");
return fCutResult;
}
}
//_________________________________________________________________________________________________
-Bool_t AliRsnCutPID::IsSelected(TObject *obj1, TObject* /*obj2*/)
+Int_t AliRsnCutPID::RealisticPID(AliRsnDaughter * const daughter, Double_t &prob)
{
//
-// Cut checker.
+// Combines the weights collected from chosen detectors with the priors
+// and gives the corresponding particle with the largest probability,
+// in terms of the AliPID particle type enumeration.
+// The argument, passed by reference, gives the corresponding probability,
+// since the cut could require that it is larger than a threshold.
//
- // coherence check
- if (!AliRsnCut::TargetOK(obj1))
+ // try to compute the weights
+ if (!ComputeWeights(daughter))
{
- AliError(Form("Wrong target. Skipping cut", GetName()));
- return kTRUE;
+ prob = -1.0;
+ return AliPID::kUnknown;
}
- // convert the object into the unique correct type
- AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(obj1);
-
- // try to compute the weights
- if (!ComputeWeights(daughter)) return kFALSE;
-
- // combine with priors and get the majority
+ // combine with priors and normalize
Int_t i;
Double_t sum = 0.0, w[AliPID::kSPECIES];
for (i = 0; i < AliPID::kSPECIES; i++)
if (sum <= 0.0)
{
AliError("Sum = 0");
- return kFALSE;
+ prob = -1.0;
+ return AliPID::kUnknown;
}
for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
// find the largest probability and related PID
- // and assign them to the mother class members which
- // are checked for the cut
- fCutValueI = 0;
- fCutValueD = w[0];
+ Int_t ibest = 0;
+ prob = w[0];
for (i = 1; i < AliPID::kSPECIES; i++)
{
- if (w[i] > fCutValueD)
+ if (w[i] > prob)
{
- fCutValueD = w[i];
- fCutValueI = i;
+ prob = w[i];
+ ibest = i;
}
}
- // if the best probability is too small, the cut is failed anyway
- if (!OkRangeD()) return kFALSE;
+ // return the value, while the probability
+ // will be consequentially set
+ return ibest;
+}
+
+//_________________________________________________________________________________________________
+Int_t AliRsnCutPID::PerfectPID(AliRsnDaughter * const daughter)
+{
+//
+// If MC is present, retrieve the particle corresponding to the used track
+// (using the fRefMC data member) and then return the true particle type
+// taken from the PDG code of the reference particle itself, converted
+// into the enumeration from AliPID object.
+//
+
+ // works only if the MC is present
+ if (!daughter->GetRefMC()) return AliPID::kUnknown;
+
+ // get the PDG code of the particle
+ TParticle *part = daughter->GetRefMC()->Particle();
+ Int_t pdg = part->GetPdgCode();
+
+ // loop over all species listed in AliPID to find the match
+ Int_t i;
+ for (i = 0; i < AliPID::kSPECIES; i++)
+ {
+ if (AliPID::ParticleCode(i) == pdg) return i;
+ }
+
+ return AliPID::kUnknown;
+}
+
+//_________________________________________________________________________________________________
+Bool_t AliRsnCutPID::IsSelected(TObject *obj1, TObject* /*obj2*/)
+{
+//
+// Cut checker.
+//
+
+ // coherence check
+ if (!AliRsnCut::TargetOK(obj1))
+ {
+ AliError(Form("Wrong target. Skipping cut", GetName()));
+ return kTRUE;
+ }
+
+ // convert the object into the unique correct type
+ AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(obj1);
- // if the best probability is OK, the cut is passed
- // if it correspond to the right particle
- return OkValue();
+ // depending on the PID type, recalls the appropriate method:
+ // in case of perfect PID, checks only if the PID type is
+ // corresponding to the request in the cut (fMinI)
+ // while in case of realistic PID checks also the probability
+ // to be within the required interval
+ if (fPerfect)
+ {
+ fCutValueI = PerfectPID(daughter);
+ return OkValueI();
+ }
+ else
+ {
+ fCutValueI = RealisticPID(daughter, fCutValueD);
+ return OkValueI() && OkRangeD();
+ }
}
+//__________________________________________________________________________________________________
void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
{
//
void IncludeDetector(EDetector det, Double_t threshold = 0., Bool_t goAbove = kTRUE);
void ExcludeDetector(EDetector det) {if (CheckBounds(det)) fUseDetector[det] = kFALSE;}
+
+ Bool_t ComputeWeights(AliRsnDaughter *daughter);
+ Int_t RealisticPID(AliRsnDaughter * const daughter, Double_t &prob);
+ Int_t PerfectPID(AliRsnDaughter * const daughter);
+ Double_t GetWeight(Int_t i) {if (i>=0&&i<AliPID::kSPECIES) return fWeight[i]; return 0.0;}
virtual Bool_t IsSelected(TObject *obj1, TObject *obj2 = 0x0);
Bool_t CheckBounds(EDetector det) const {return (det >= kITS && det < kDetectors);}
Bool_t CheckThreshold(EDetector det, Double_t value);
- Bool_t ComputeWeights(AliRsnDaughter *daughter);
Double_t fPrior[AliPID::kSPECIES]; // prior probability
Double_t fWeight[AliPID::kSPECIES]; // PID weights used for combinations
return kTRUE;
}
+//_____________________________________________________________________________
+Bool_t AliRsnEvent::SetDaughterMC(AliRsnDaughter &out, Int_t i)
+{
+//
+// Using the second argument, retrieves the i-th object in the
+// appropriate sample and sets the firs reference object
+// in order to point to that.
+// If a MonteCarlo information is provided, sets the useful informations from there,
+// and in case of a V0, sets the 'label' data member only when the two daughters
+// of the V0 point to the same mother.
+// Returns kFALSE whenever the operation fails (out of range, NULL references).
+//
+
+ if (!fRefMC)
+ {
+ out.SetBad();
+ return kFALSE;
+ }
+
+ if (i >= fRefMC->GetNumberOfTracks())
+ {
+ out.SetBad();
+ return kFALSE;
+ }
+
+ AliMCParticle *track = (AliMCParticle*)fRef->GetTrack(i);
+ if (!track)
+ {
+ out.SetBad();
+ return kFALSE;
+ }
+ else
+ {
+ out.SetRef(track);
+ out.SetRefMC(track);
+ out.SetLabel(i);
+ out.SetGood();
+ }
+
+ AliStack *stack = fRefMC->Stack();
+ TParticle *part = track->Particle();
+ if (part)
+ {
+ Int_t imum = part->GetFirstMother();
+ if (imum >= 0 && imum <= stack->GetNtrack())
+ {
+ TParticle *mum = stack->Particle(imum);
+ if (mum) out.SetMotherPDG(TMath::Abs(mum->GetPdgCode()));
+ }
+ }
+
+ return kTRUE;
+}
+
//_____________________________________________________________________________
AliRsnDaughter AliRsnEvent::GetDaughter(Int_t i, AliRsnDaughter::ERefType type)
{
return out;
}
+//_____________________________________________________________________________
+AliRsnDaughter AliRsnEvent::GetDaughterMC(Int_t i)
+{
+//
+// Return an AliRsnDaughter taken from this event,
+// with all additional data members well set.
+//
+
+ AliRsnDaughter out;
+ SetDaughterMC(out, i);
+
+ return out;
+}
+
//_____________________________________________________________________________
Int_t AliRsnEvent::GetMultiplicity()
{
Int_t GetMultiplicity();
Bool_t SetDaughter(AliRsnDaughter &daughter, Int_t index, AliRsnDaughter::ERefType type = AliRsnDaughter::kTrack);
+ Bool_t SetDaughterMC(AliRsnDaughter &daughter, Int_t index);
AliRsnDaughter GetDaughter(Int_t i, AliRsnDaughter::ERefType type = AliRsnDaughter::kTrack);
+ AliRsnDaughter GetDaughterMC(Int_t i);
Int_t SelectLeadingParticle(Double_t ptMin = 0.0, AliRsnCutPID *cutPID = 0x0);
Int_t GetLeadingParticleID() {return fLeading;}
//________________________________________________________________________________________
AliRsnFunction::AliRsnFunction(Bool_t useTH1) :
- TNamed(),
+ TObject(),
fPairDef(0x0),
fAxisList("AliRsnValue", 0),
fPair(0x0),
//________________________________________________________________________________________
AliRsnFunction::AliRsnFunction(const AliRsnFunction ©) :
- TNamed(copy),
+ TObject(copy),
fPairDef(copy.fPairDef),
fAxisList(copy.fAxisList),
fPair(copy.fPair),
// Assignment operator.
//
- SetName(copy.GetName());
- SetTitle(copy.GetTitle());
+ //SetName(copy.GetName());
+ //SetTitle(copy.GetTitle());
fPairDef = copy.fPairDef;
fPair = copy.fPair;
class AliRsnMother;
class AliRsnPairDef;
-class AliRsnFunction : public TNamed
+class AliRsnFunction : public TObject
{
public:
}
//_____________________________________________________________________________
-Double_t AliRsnMother::ThetaStar(Bool_t first, Bool_t useMC)
+Double_t AliRsnMother::CosThetaStar(Bool_t first, Bool_t useMC)
{
-//
-// Returns the theta* as the angle of the first daughter
-// w.r. to the mother momentum, in its rest frame
-//
+ TLorentzVector mother = (useMC ? fSumMC : fSum);
+ TLorentzVector daughter0 = (first ? fDaughter[0]->P() : fDaughter[1]->P());
+ TLorentzVector daughter1 = (first ? fDaughter[1]->P() : fDaughter[0]->P());
+ TVector3 momentumM (mother.Vect());
+ //cout << mother.Vect().X() << ' ' << mother.Vect().Y() << ' ' << mother.Vect().Z() << endl;
+ TVector3 normal (mother.Y()/momentumM.Mag(), -mother.X()/momentumM.Mag(), 0.0);
+ //cout << momentumM.Mag() << ' ' << normal.X() << ' ' << normal.Y() << ' ' << normal.Z() << endl;
+
+// Computes first the invariant mass of the mother
+
+ Double_t mass0 = fDaughter[0]->P().M();
+ Double_t mass1 = fDaughter[1]->P().M();
+ Double_t p0 = daughter0.Vect().Mag();
+ Double_t p1 = daughter1.Vect().Mag();
+ Double_t E0 = TMath::Sqrt(mass0*mass0+p0*p0);
+ Double_t E1 = TMath::Sqrt(mass1*mass1+p1*p1);
+ Double_t MotherMass = TMath::Sqrt((E0+E1)*(E0+E1)-(p0*p0+2.0*daughter0.Vect().Dot(daughter1.Vect())+p1*p1));
+ MotherMass = fSum.M();
+
+// Computes components of beta
+
+ Double_t betaX = -mother.X()/mother.E(); //TMath::Sqrt(momentumM.Mag()*momentumM.Mag()+MotherMass*MotherMass);
+ Double_t betaY = -mother.Y()/mother.E(); //TMath::Sqrt(momentumM.Mag()*momentumM.Mag()+MotherMass*MotherMass);
+ Double_t betaZ = -mother.Z()/mother.E(); //TMath::Sqrt(momentumM.Mag()*momentumM.Mag()+MotherMass*MotherMass);
+
+// Computes Lorentz transformation of the momentum of the first daughter
+// into the rest frame of the mother and theta*
+ //cout << "Beta = " << betaX << ' ' << betaY << ' ' << betaZ << endl;
+
+ daughter0.Boost(betaX,betaY,betaZ);
+ TVector3 momentumD = daughter0.Vect();
+
+ Double_t cosThetaStar = normal.Dot(momentumD)/momentumD.Mag();
- TLorentzVector &mother = (useMC ? fSumMC : fSum);
- TLorentzVector &daughter = (first ? fDaughter[0]->P() : fDaughter[1]->P());
-
- Double_t theta = mother.Vect().Theta();
- Double_t phi = mother.Vect().Phi();
- Double_t beta = mother.Beta();
-
- Double_t betax = beta * TMath::Sin(theta) * TMath::Cos(phi);
- Double_t betay = beta * TMath::Sin(theta) * TMath::Sin(phi);
- Double_t betaz = beta * TMath::Cos(theta);
- Double_t gamx = 1.0 / TMath::Sqrt(1.0 - betax*betax);
- Double_t gamy = 1.0 / TMath::Sqrt(1.0 - betay*betay);
- Double_t gamz = 1.0 / TMath::Sqrt(1.0 - betaz*betaz);
-
- Double_t pstarx = gamx * (daughter.Vect().X() - betax * daughter.E());
- Double_t pstary = gamy * (daughter.Vect().Y() - betay * daughter.E());
- Double_t pstarz = gamz * (daughter.Vect().Z() - betaz * daughter.E());
-
- TVector3 dstar(pstarx, pstary, pstarz);
- TVector3 vnorm(mother.Vect().Y() / mother.Perp(), mother.Vect().X() / mother.Perp(), 0.0);
-
- Double_t cosThetaStar = dstar.Dot(vnorm) / dstar.Mag();
-
return cosThetaStar;
+
}
//_____________________________________________________________________________
TLorentzVector& RefMC() {return fRefMC;}
Double_t OpeningAngle(Bool_t mc = kFALSE) const {if (fDaughter[0] && fDaughter[1]) return fDaughter[0]->P(mc).Angle(fDaughter[1]->P(mc).Vect()); return 1E6;}
Double_t AngleTo(AliRsnDaughter track, Bool_t mc = kFALSE) const {return fSum.Angle(track.P(mc).Vect());}
- Double_t ThetaStar(Bool_t first = kTRUE, Bool_t useMC = kFALSE);
+ Double_t CosThetaStar(Bool_t first = kTRUE, Bool_t useMC = kFALSE);
AliRsnDaughter* GetDaughter(const Int_t &index) const {if (index==0||index==1) return fDaughter[index]; return 0x0;}
}
//_____________________________________________________________________________
-Bool_t AliRsnValue::Eval(AliRsnMother *mother, const AliRsnPairDef *pairDef, AliRsnEvent *event)
+Bool_t AliRsnValue::Eval(AliRsnMother * const mother, AliRsnPairDef * const pairDef, AliRsnEvent * const event)
{
//
// Evaluation of the required value.
mother->SetDefaultMass(mass);
fValue = mother->Ref().Rapidity();
break;
+ case kPairCosThetaStar:
+ fValue = mother->CosThetaStar();
+ break;
case kPairCosThetaStar1:
- fValue = TMath::Cos(mother->ThetaStar(kTRUE, kFALSE));
+ //fValue = TMath::Cos(mother->ThetaStar(kTRUE, kFALSE));
break;
case kPairCosThetaStar2:
- fValue = TMath::Cos(mother->ThetaStar(kFALSE, kFALSE));
+ //fValue = TMath::Cos(mother->ThetaStar(kFALSE, kFALSE));
break;
case kPairCosThetaStarMC1:
- fValue = TMath::Cos(mother->ThetaStar(kTRUE, kTRUE));
+ //fValue = TMath::Cos(mother->ThetaStar(kTRUE, kTRUE));
break;
case kPairCosThetaStarMC2:
- fValue = TMath::Cos(mother->ThetaStar(kFALSE, kTRUE));
+ //fValue = TMath::Cos(mother->ThetaStar(kFALSE, kTRUE));
break;
case kEventMult:
if (!event) fValue = 0.0;
kPairEta,
kPairMt,
kPairY,
+ kPairCosThetaStar,
kPairCosThetaStar1,
kPairCosThetaStar2,
kPairCosThetaStarMC1,
void SetBins(Int_t n, Double_t min, Double_t max);
void SetBins(Double_t min, Double_t max, Double_t step);
- Bool_t Eval(AliRsnMother *mother, const AliRsnPairDef *pairDef, AliRsnEvent *event);
+ Bool_t Eval(AliRsnMother * const mother, AliRsnPairDef * const pairDef, AliRsnEvent * const event);
private:
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff