2 // *** Class AliRsnLoopPair ***
4 // "Core" method for defining the work on a pari of particles.
5 // For one analysis, one must setup one of this for each pair he wants to analyze,
6 // adding to it all analysis which he desires to do.
7 // Here he defines the cuts, and the particle types and charges, and can add
8 // functions which do different operations on the same pair, and some binning
9 // with respect to some kinematic variables (eta, momentum)
11 // authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
12 // M. Vala (email: martin.vala@cern.ch)
15 #include <Riostream.h>
17 #include <TEntryList.h>
21 #include "AliRsnEvent.h"
22 #include "AliRsnPairDef.h"
23 #include "AliRsnMother.h"
24 #include "AliRsnCutSet.h"
25 #include "AliRsnDaughterSelector.h"
27 #include "AliRsnLoopPair.h"
29 ClassImp(AliRsnLoopPair)
31 //_____________________________________________________________________________
32 AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMixed) :
33 AliRsnLoop(name, isMixed),
44 // Default constructor
51 //_____________________________________________________________________________
52 AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair ©) :
54 fTrueMC(copy.fTrueMC),
55 fOnlyTrue(copy.fOnlyTrue),
56 fUseMCRef(copy.fUseMCRef),
57 fCheckDecay(copy.fCheckDecay),
58 fRangeY(copy.fRangeY),
59 fPairDef(copy.fPairDef),
60 fPairCuts(copy.fPairCuts),
67 fListID[0] = copy.fListID[0];
68 fListID[1] = copy.fListID[1];
71 //_____________________________________________________________________________
72 AliRsnLoopPair &AliRsnLoopPair::operator=(const AliRsnLoopPair ©)
75 // Assignment operator
78 AliRsnLoop::operator=(copy);
81 fTrueMC = copy.fTrueMC;
82 fOnlyTrue = copy.fOnlyTrue;
83 fUseMCRef = copy.fUseMCRef;
84 fCheckDecay = copy.fCheckDecay;
85 fRangeY = copy.fRangeY;
86 fPairDef = copy.fPairDef;
87 fPairCuts = copy.fPairCuts;
88 fMother = copy.fMother;
89 fListID[0] = copy.fListID[0];
90 fListID[1] = copy.fListID[1];
95 //_____________________________________________________________________________
96 AliRsnLoopPair::~AliRsnLoopPair()
103 //_____________________________________________________________________________
104 void AliRsnLoopPair::Print(Option_t * /*option*/) const
107 // Prints info about pair
113 //_____________________________________________________________________________
114 Bool_t AliRsnLoopPair::Init(const char *prefix, TList *list)
117 // Initialization function.
118 // Loops on all functions and eventual the ntuple, to initialize output objects.
121 // assign data members relations
122 fMother.SetDaughter(0, &fDaughter[0]);
123 fMother.SetDaughter(1, &fDaughter[1]);
124 AliInfo(Form("[%s] Initialization", GetName()));
126 TString name(prefix);
129 // if (IsMixed()) name.Append("_mix");
131 return AliRsnLoop::Init(name.Data(), list);
134 //_____________________________________________________________________________
135 Int_t AliRsnLoopPair::DoLoop
136 (AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
140 // Computes what is needed from passed events.
141 // Returns the number of pairs successfully processed.
145 AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
146 if (!evMix || !selMix) {
147 AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
151 AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
155 fMother.SetRefEvent(evMain);
158 if (!OkEvent(evMain)) {
159 AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
162 if (!OkEvent(evMix)) {
163 AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
167 // if it is required to loop over True MC, do this here and skip the rest of the method
168 if (fTrueMC) return LoopTrueMC(evMain);
170 Int_t i0, i1, start, npairs = 0;
172 TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
173 TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
174 if (!list0 || !list1) {
175 AliError("Can't process NULL lists");
178 AliDebugClass(3, Form("[%s]: list counts: %lld, %lld", GetName(), list0->GetN(), list1->GetN()));
179 if (!list0->GetN() || !list1->GetN()) {
180 AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
184 TObjArrayIter next(&fOutputs);
185 AliRsnListOutput *out = 0x0;
186 Long64_t iEntry1,iEntry2;
187 for (i0 = 0; i0 < list0->GetN(); i0++) {
188 iEntry1 = list0->GetEntry(i0);
189 evMain->SetDaughter(fDaughter[0], iEntry1,fUseMCRef);
190 fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
192 if (!fIsMixed && list0 == list1) start = i0 + 1;
193 for (i1 = start; i1 < list1->GetN(); i1++) {
194 iEntry2 = list1->GetEntry(i1);
195 if (iEntry1 == iEntry2) continue;
196 AliDebugClass(4, Form("Checking entries pair: %d (%lld) with %d (%lld)", i0, iEntry1, i1, iEntry2));
197 evMix->SetDaughter(fDaughter[1], iEntry2,fUseMCRef);
198 fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
199 fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
200 fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
201 fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
202 fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
203 // check rapidity range
204 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
205 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
210 if (!IsTrueMother()) {
211 AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
217 if (!fPairCuts->IsSelected(&fMother)) {
218 AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
224 while ( (out = (AliRsnListOutput *)next()) ) {
225 if (out->Fill(&fMother)) npairs++;
226 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
234 //_____________________________________________________________________________
235 Bool_t AliRsnLoopPair::IsTrueMother()
238 // Checks to see if the mother comes from a true resonance.
239 // It is triggered by the 'SetOnlyTrue()' function
243 // daughters have same mother with the right PDG code
244 Int_t commonPDG = fMother.CommonMother();
245 if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
246 AliDebugClass(1, "Found a true mother");
249 // checks if daughter have the right particle type
250 // (activated by fCheckDecay)
252 AliRsnDaughterDef &def1 = fPairDef->GetDef1();
253 AliRsnDaughterDef &def2 = fPairDef->GetDef2();
254 if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
255 if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
257 AliDebugClass(1, "Decay products match");
262 //__________________________________________________________________________________________________
263 Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
266 // Calls the appropriate assignment method
270 fMother.SetDaughter(0, &fDaughter[0]);
271 fMother.SetDaughter(1, &fDaughter[1]);
272 fMother.SetRefEvent(rsnEvent);
273 fDaughter[0].SetOwnerEvent(rsnEvent);
274 fDaughter[1].SetOwnerEvent(rsnEvent);
276 if (rsnEvent->IsESD())
277 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
278 else if (rsnEvent->IsAOD())
279 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
281 AliError("Unrecognized input event");
286 //__________________________________________________________________________________________________
287 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
290 // Gets a particle in the MC event and try to assign it to the mother.
291 // If it has two daughters, retrieve them and assign also them.
292 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
293 // for ESD and AOD, if available
295 // Implementation for ESD inputs
298 AliMCEvent *mc = rsnEvent->GetRefMCESD();
299 AliStack *stack = mc->Stack();
300 AliMCParticle *mother = (AliMCParticle *)mc->GetTrack(ipart);
301 TParticle *motherP = mother->Particle();
302 Int_t ntracks = stack->GetNtrack();
304 // check PDG code and exit if it is wrong
305 if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
307 // check number of daughters and exit if it is not 2
308 if (motherP->GetNDaughters() < 2) return kFALSE;
309 //if (!stack->IsPhysicalPrimary(ipart)) return kFALSE;
312 // check distance from primary vertex
313 TLorentzVector vprod;
314 motherP->ProductionVertex(vprod);
315 if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
316 AliDebugClass(1, "Distant production vertex");
321 // get the daughters and check their PDG code and charge:
322 // if they match one of the pair daughter definitions,
323 // assign them as MC reference of the 'fDaughter' objects
324 fDaughter[0].Reset();
325 fDaughter[1].Reset();
326 Int_t index[2] = {motherP->GetDaughter(0), motherP->GetDaughter(1)};
329 AliMCParticle *daughter = 0x0;
330 for (i = 0; i < 2; i++) {
331 // check index for stack
332 if (index[i] < 0 || index[i] > ntracks) {
333 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
336 // get daughter and its PDG and charge
337 daughter = (AliMCParticle *)mc->GetTrack(index[i]);
338 pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
339 charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
340 // check if it matches one definition
341 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
342 fDaughter[0].SetGood();
343 fDaughter[0].SetRefMC(daughter);
344 fDaughter[0].SetLabel(index[i]);
345 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
346 fDaughter[1].SetGood();
347 fDaughter[1].SetRefMC(daughter);
348 fDaughter[1].SetLabel(index[i]);
352 // return success if both daughters were assigned
353 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
356 fDaughter[0].Reset();
357 fDaughter[1].Reset();
362 //__________________________________________________________________________________________________
363 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
366 // Gets a particle in the MC event and try to assign it to the mother.
367 // If it has two daughters, retrieve them and assign also them.
368 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
369 // for ESD and AOD, if available
371 // Implementation for AOD inputs
374 AliAODEvent *aod = rsnEvent->GetRefAOD();
375 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
376 AliAODMCParticle *mother = (AliAODMCParticle *)listAOD->At(ipart);
378 // check PDG code and exit if it is wrong
379 if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
381 // check number of daughters and exit if it is not 2
382 if (mother->GetNDaughters() < 2) return kFALSE;
383 if (!mother->IsPrimary()) return kFALSE;
386 // check distance from primary vertex
387 Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
388 Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
389 if (dv > fMaxDistPV) {
390 AliDebugClass(1, "Distant production vertex");
395 // get the daughters and check their PDG code and charge:
396 // if they match one of the pair daughter definitions,
397 // assign them as MC reference of the 'fDaughter' objects
398 fDaughter[0].Reset();
399 fDaughter[1].Reset();
400 Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
401 Int_t ntracks = listAOD->GetEntriesFast();
404 AliAODMCParticle *daughter = 0x0;
405 for (i = 0; i < 2; i++) {
406 // check index for stack
407 if (index[i] < 0 || index[i] > ntracks) {
408 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
411 // get daughter and its PDG and charge
412 daughter = (AliAODMCParticle *)listAOD->At(index[i]);
413 pdg = TMath::Abs(daughter->GetPdgCode());
414 charge = (Short_t)(daughter->Charge() / 3);
415 // check if it matches one definition
416 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
417 fDaughter[0].SetGood();
418 fDaughter[0].SetRefMC(daughter);
419 fDaughter[0].SetLabel(index[i]);
420 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
421 fDaughter[1].SetGood();
422 fDaughter[1].SetRefMC(daughter);
423 fDaughter[1].SetLabel(index[i]);
427 // return success if both daughters were assigned
428 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
431 fDaughter[0].Reset();
432 fDaughter[1].Reset();
437 //_____________________________________________________________________________
438 Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
441 // Loop on event and fill containers
444 // check presence of MC reference
445 if (!rsn->GetRefMC()) {
446 AliError("Need a MC to compute efficiency");
451 // must be ESD or AOD, and then use a bool to know in the rest
452 if (!rsn->IsESD() && !rsn->IsAOD()) {
453 AliError("Need to process ESD or AOD input");
457 // retrieve the MC primary vertex position
458 // and do some additional coherence checks
459 //Double_t fVertex[3] = {0.0, 0.0, 0.0};
462 //TArrayF mcVertex(3);
463 //AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
464 //genEH->PrimaryVertex(mcVertex);
465 //for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
466 npart = rsn->GetRefMCESD()->GetNumberOfTracks();
468 //for (i = 0; i < 3; i++) fVertex[i] = 0.0;
469 AliAODEvent *aod = rsn->GetRefMCAOD();
470 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
471 if (listAOD) npart = listAOD->GetEntries();
472 //AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
473 //if (mcH) mcH->GetVertex(fVertex);
476 // check number of particles
478 AliInfo("Empty event");
483 Int_t ipart, count = 0;
484 AliRsnDaughter check;
486 TObjArrayIter next(&fOutputs);
487 AliRsnListOutput *out = 0x0;
489 // loop over particles
490 for (ipart = 0; ipart < npart; ipart++) {
491 // check i-th particle
492 if (!AssignMotherAndDaughters(rsn, ipart)) continue;
493 // if assignment was successful, for first step, use MC info
494 fDaughter[0].SetRef(fDaughter[0].GetRefMC());
495 fDaughter[1].SetRef(fDaughter[1].GetRefMC());
496 fMother.SetRefEvent(rsn);
497 fMother.ComputeSum(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
498 // check rapidity range
499 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
500 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
505 while ( (out = (AliRsnListOutput *)next()) ) {
506 if (out->Fill(&fMother)) count++;
507 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));