2 // Class AliRsnLoopEffPair
4 // Inherits from basic AliRsnLoopEff for efficiency,
5 // and computed efficiencies for single-tracks
7 // author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
10 #include <Riostream.h>
13 #include "AliGenEventHeader.h"
14 #include "AliAODMCHeader.h"
16 #include "AliRsnEvent.h"
17 #include "AliRsnCutManager.h"
18 #include "AliRsnDaughterDef.h"
19 #include "AliRsnPairDef.h"
20 #include "AliRsnLoopEffPair.h"
22 ClassImp(AliRsnLoopEffPair)
24 //_____________________________________________________________________________
25 AliRsnLoopEffPair::AliRsnLoopEffPair(const char *name, AliRsnPairDef *def) :
26 AliRsnLoopEff(name, 1),
31 // Default constructor.
32 // Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
36 //_____________________________________________________________________________
37 AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair& copy) :
47 //_____________________________________________________________________________
48 AliRsnLoopEffPair& AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair& copy)
51 // Assignment operator.
52 // Owned data members are meaningless for this operator.
55 AliRsnLoopEff::operator=(copy);
63 //__________________________________________________________________________________________________
64 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
67 // Calls the appropriate assignment method
71 fMother.SetDaughter(0, &fDaughter[0]);
72 fMother.SetDaughter(1, &fDaughter[1]);
73 fMother.SetRefEvent(rsnEvent);
74 fDaughter[0].SetOwnerEvent(rsnEvent);
75 fDaughter[1].SetOwnerEvent(rsnEvent);
77 if (rsnEvent->IsESD())
78 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
79 else if (rsnEvent->IsAOD())
80 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
82 AliError("Unrecognized input event");
87 //__________________________________________________________________________________________________
88 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
91 // Gets a particle in the MC event and try to assign it to the mother.
92 // If it has two daughters, retrieve them and assign also them.
93 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
94 // for ESD and AOD, if available
96 // Implementation for ESD inputs
99 AliMCEvent *mc = rsnEvent->GetRefMCESD();
100 AliStack *stack = mc->Stack();
101 AliMCParticle *mother = (AliMCParticle*)mc->GetTrack(ipart);
102 TParticle *motherP = mother->Particle();
103 Int_t ntracks = stack->GetNtrack();
105 // check PDG code and exit if it is wrong
106 if (TMath::Abs(motherP->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
108 // check number of daughters and exit if it is not 2
109 if (motherP->GetNDaughters() < 2) return kFALSE;
111 // check distance from primary vertex
112 TLorentzVector vprod;
113 motherP->ProductionVertex(vprod);
114 if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
115 AliDebugClass(1, "Distant production vertex");
119 // get the daughters and check their PDG code and charge:
120 // if they match one of the pair daughter definitions,
121 // assign them as MC reference of the 'fDaughter' objects
122 fDaughter[0].Reset();
123 fDaughter[1].Reset();
124 Int_t index[2] = {motherP->GetFirstDaughter(), motherP->GetLastDaughter()};
127 AliMCParticle *daughter = 0x0;
128 for (i = 0; i < 2; i++) {
129 // check index for stack
130 if (index[i] < 0 || index[i] > ntracks) {
131 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
134 // get daughter and its PDG and charge
135 daughter = (AliMCParticle*)mc->GetTrack(index[i]);
136 pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
137 charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
138 // check if it matches one definition
139 if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
140 fDaughter[0].SetGood();
141 fDaughter[0].SetRefMC(daughter);
142 fDaughter[0].SetLabel(index[i]);
143 } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
144 fDaughter[1].SetGood();
145 fDaughter[1].SetRefMC(daughter);
146 fDaughter[1].SetLabel(index[i]);
150 // return success if both daughters were assigned
151 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
154 fDaughter[0].Reset();
155 fDaughter[1].Reset();
160 //__________________________________________________________________________________________________
161 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
164 // Gets a particle in the MC event and try to assign it to the mother.
165 // If it has two daughters, retrieve them and assign also them.
166 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
167 // for ESD and AOD, if available
169 // Implementation for AOD inputs
172 AliAODEvent *aod = rsnEvent->GetRefAOD();
173 TClonesArray *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
174 AliAODMCParticle *mother = (AliAODMCParticle*)listAOD->At(ipart);
175 Int_t ntracks = listAOD->GetEntries();
177 // check PDG code and exit if it is wrong
178 if (TMath::Abs(mother->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
180 // check number of daughters and exit if it is not 2
181 if (mother->GetNDaughters() < 2) return kFALSE;
183 // check distance from primary vertex
184 Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
185 Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
186 if (dv > fMaxDistPV) {
187 AliDebugClass(1, "Distant production vertex");
191 // get the daughters and check their PDG code and charge:
192 // if they match one of the pair daughter definitions,
193 // assign them as MC reference of the 'fDaughter' objects
194 fDaughter[0].Reset();
195 fDaughter[1].Reset();
196 Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
199 AliAODMCParticle *daughter = 0x0;
200 for (i = 0; i < 2; i++) {
201 // check index for stack
202 if (index[i] < 0 || index[i] > ntracks) {
203 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
206 // get daughter and its PDG and charge
207 daughter = (AliAODMCParticle*)listAOD->At(index[i]);
208 pdg = TMath::Abs(daughter->GetPdgCode());
209 charge = (Short_t)(daughter->Charge() / 3);
210 // check if it matches one definition
211 if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
212 fDaughter[0].SetGood();
213 fDaughter[0].SetRefMC(daughter);
214 fDaughter[0].SetLabel(index[i]);
215 } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
216 fDaughter[1].SetGood();
217 fDaughter[1].SetRefMC(daughter);
218 fDaughter[1].SetLabel(index[i]);
222 // return success if both daughters were assigned
223 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
226 fDaughter[0].Reset();
227 fDaughter[1].Reset();
232 //_____________________________________________________________________________
233 Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent *rsn, AliRsnDaughterSelector*, AliRsnEvent*, AliRsnDaughterSelector*)
236 // Loop on event and fill containers
240 if (!OkEvent(rsn)) return 0;
243 fOutput = (AliRsnListOutput*)fOutputs[0];
245 // check presence of MC reference
246 if (!rsn->GetRefMC()) {
247 AliError("Need a MC to compute efficiency");
251 // check presence of event
252 if (!rsn->GetRef()) {
253 AliError("Need an event to compute efficiency");
258 // must be ESD or AOD, and then use a bool to know in the rest
259 if (!rsn->IsESD() && !rsn->IsAOD()) {
260 AliError("Need to process ESD or AOD input");
264 // retrieve the MC primary vertex position
265 // and do some additional coherence checks
269 AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
270 genEH->PrimaryVertex(mcVertex);
271 for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
272 npart = rsn->GetRefMCESD()->GetNumberOfTracks();
274 for (i = 0; i < 3; i++) fVertex[i] = 0.0;
275 AliAODEvent *aod = rsn->GetRefMCAOD();
276 TClonesArray *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
277 if (listAOD) npart = listAOD->GetEntries();
278 AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
279 if (mcH) mcH->GetVertex(fVertex);
282 // check number of particles
284 AliInfo("Empty event");
289 Int_t ipart, istep, count = 0, nsteps = fSteps.GetEntries();
290 Int_t ntracks = rsn->GetAbsoluteSum();
291 AliRsnDaughter check;
293 // loop over particles
294 for (ipart = 0; ipart < npart; ipart++) {
295 // check i-th particle
296 if (!AssignMotherAndDaughters(rsn, ipart)) continue;
297 // if assignment was successful, for first step, use MC info
298 fDaughter[0].SetRef(fDaughter[0].GetRefMC());
299 fDaughter[1].SetRef(fDaughter[1].GetRefMC());
300 fMother.SetRefEvent(rsn);
301 fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
302 fOutput->Fill(&fMother, 0);
304 // for each further step, try to find two tracks which pass the related cuts
305 for (istep = 0; istep < nsteps; istep++) {
306 AliRsnCutManager *cuts = (AliRsnCutManager*)fSteps[istep];
307 fDaughter[0].SetBad();
308 fDaughter[1].SetBad();
309 for (i = 0; i < ntracks; i++) {
310 rsn->SetDaughter(check, i);
311 if (!cuts->PassCommonDaughterCuts(&check)) continue;
312 if (TMath::Abs(check.GetLabel()) == fDaughter[0].GetLabel()) {
313 if (!cuts->PassDaughter1Cuts(&check)) continue;
314 fDaughter[0].SetRef(check.GetRef());
315 fDaughter[0].SetGood();
316 } else if (TMath::Abs(check.GetLabel()) == fDaughter[1].GetLabel()) {
317 if (!cuts->PassDaughter2Cuts(&check)) continue;
318 fDaughter[1].SetRef(check.GetRef());
319 fDaughter[1].SetGood();
321 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
322 fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
323 if (cuts->PassMotherCuts(&fMother)) {
324 fOutput->Fill(&fMother, 1 + istep);