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.Prepend("mix_");
130 if (IsMixed()) name.Append("_mix");
132 return AliRsnLoop::Init(name.Data(), list);
135 //_____________________________________________________________________________
136 Int_t AliRsnLoopPair::DoLoop
137 (AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
141 // Computes what is needed from passed events.
142 // Returns the number of pairs successfully processed.
146 AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
147 if (!evMix || !selMix) {
148 AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
152 AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
156 fMother.SetRefEvent(evMain);
159 if (!OkEvent(evMain)) {
160 AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
163 if (!OkEvent(evMix)) {
164 AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
168 // if it is required to loop over True MC, do this here and skip the rest of the method
169 if (fTrueMC) return LoopTrueMC(evMain);
171 Int_t i0, i1, start, npairs = 0;
173 TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
174 TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
175 if (!list0 || !list1) {
176 AliError("Can't process NULL lists");
179 AliDebugClass(3, Form("[%s]: list counts: %lld, %lld", GetName(), list0->GetN(), list1->GetN()));
180 if (!list0->GetN() || !list1->GetN()) {
181 AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
185 TObjArrayIter next(&fOutputs);
186 AliRsnListOutput *out = 0x0;
187 Long64_t iEntry1,iEntry2;
188 for (i0 = 0; i0 < list0->GetN(); i0++) {
189 iEntry1 = list0->GetEntry(i0);
190 evMain->SetDaughter(fDaughter[0], iEntry1,fUseMCRef);
191 fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
193 if (!fIsMixed && list0 == list1) start = i0 + 1;
194 for (i1 = start; i1 < list1->GetN(); i1++) {
195 iEntry2 = list1->GetEntry(i1);
196 if (iEntry1 == iEntry2) continue;
197 AliDebugClass(4, Form("Checking entries pair: %d (%lld) with %d (%lld)", i0, iEntry1, i1, iEntry2));
198 evMix->SetDaughter(fDaughter[1], iEntry2,fUseMCRef);
199 fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
200 fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
201 fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
202 fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
203 fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
204 // check rapidity range
205 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
206 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
211 if (!IsTrueMother()) {
212 AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
218 if (!fPairCuts->IsSelected(&fMother)) {
219 AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
225 while ( (out = (AliRsnListOutput *)next()) ) {
226 if (out->Fill(&fMother)) npairs++;
227 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
235 //_____________________________________________________________________________
236 Bool_t AliRsnLoopPair::IsTrueMother()
239 // Checks to see if the mother comes from a true resonance.
240 // It is triggered by the 'SetOnlyTrue()' function
244 // daughters have same mother with the right PDG code
245 Int_t commonPDG = fMother.CommonMother();
246 if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
247 AliDebugClass(1, "Found a true mother");
250 // checks if daughter have the right particle type
251 // (activated by fCheckDecay)
253 AliRsnDaughterDef &def1 = fPairDef->GetDef1();
254 AliRsnDaughterDef &def2 = fPairDef->GetDef2();
255 if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
256 if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
258 AliDebugClass(1, "Decay products match");
263 //__________________________________________________________________________________________________
264 Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
267 // Calls the appropriate assignment method
271 fMother.SetDaughter(0, &fDaughter[0]);
272 fMother.SetDaughter(1, &fDaughter[1]);
273 fMother.SetRefEvent(rsnEvent);
274 fDaughter[0].SetOwnerEvent(rsnEvent);
275 fDaughter[1].SetOwnerEvent(rsnEvent);
277 if (rsnEvent->IsESD())
278 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
279 else if (rsnEvent->IsAOD())
280 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
282 AliError("Unrecognized input event");
287 //__________________________________________________________________________________________________
288 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
291 // Gets a particle in the MC event and try to assign it to the mother.
292 // If it has two daughters, retrieve them and assign also them.
293 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
294 // for ESD and AOD, if available
296 // Implementation for ESD inputs
299 AliMCEvent *mc = rsnEvent->GetRefMCESD();
300 AliStack *stack = mc->Stack();
301 AliMCParticle *mother = (AliMCParticle *)mc->GetTrack(ipart);
302 TParticle *motherP = mother->Particle();
303 Int_t ntracks = stack->GetNtrack();
305 // check PDG code and exit if it is wrong
306 if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
308 // check number of daughters and exit if it is not 2
309 if (motherP->GetNDaughters() < 2) return kFALSE;
310 //if (!stack->IsPhysicalPrimary(ipart)) return kFALSE;
313 // check distance from primary vertex
314 TLorentzVector vprod;
315 motherP->ProductionVertex(vprod);
316 if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
317 AliDebugClass(1, "Distant production vertex");
322 // get the daughters and check their PDG code and charge:
323 // if they match one of the pair daughter definitions,
324 // assign them as MC reference of the 'fDaughter' objects
325 fDaughter[0].Reset();
326 fDaughter[1].Reset();
327 Int_t index[2] = {motherP->GetDaughter(0), motherP->GetDaughter(1)};
330 AliMCParticle *daughter = 0x0;
331 for (i = 0; i < 2; i++) {
332 // check index for stack
333 if (index[i] < 0 || index[i] > ntracks) {
334 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
337 // get daughter and its PDG and charge
338 daughter = (AliMCParticle *)mc->GetTrack(index[i]);
339 pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
340 charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
341 // check if it matches one definition
342 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
343 fDaughter[0].SetGood();
344 fDaughter[0].SetRefMC(daughter);
345 fDaughter[0].SetLabel(index[i]);
346 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
347 fDaughter[1].SetGood();
348 fDaughter[1].SetRefMC(daughter);
349 fDaughter[1].SetLabel(index[i]);
353 // return success if both daughters were assigned
354 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
357 fDaughter[0].Reset();
358 fDaughter[1].Reset();
363 //__________________________________________________________________________________________________
364 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
367 // Gets a particle in the MC event and try to assign it to the mother.
368 // If it has two daughters, retrieve them and assign also them.
369 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
370 // for ESD and AOD, if available
372 // Implementation for AOD inputs
375 AliAODEvent *aod = rsnEvent->GetRefAOD();
376 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
377 AliAODMCParticle *mother = (AliAODMCParticle *)listAOD->At(ipart);
379 // check PDG code and exit if it is wrong
380 if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
382 // check number of daughters and exit if it is not 2
383 if (mother->GetNDaughters() < 2) return kFALSE;
384 if (!mother->IsPrimary()) return kFALSE;
387 // check distance from primary vertex
388 Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
389 Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
390 if (dv > fMaxDistPV) {
391 AliDebugClass(1, "Distant production vertex");
396 // get the daughters and check their PDG code and charge:
397 // if they match one of the pair daughter definitions,
398 // assign them as MC reference of the 'fDaughter' objects
399 fDaughter[0].Reset();
400 fDaughter[1].Reset();
401 Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
402 Int_t ntracks = listAOD->GetEntriesFast();
405 AliAODMCParticle *daughter = 0x0;
406 for (i = 0; i < 2; i++) {
407 // check index for stack
408 if (index[i] < 0 || index[i] > ntracks) {
409 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
412 // get daughter and its PDG and charge
413 daughter = (AliAODMCParticle *)listAOD->At(index[i]);
414 pdg = TMath::Abs(daughter->GetPdgCode());
415 charge = (Short_t)(daughter->Charge() / 3);
416 // check if it matches one definition
417 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
418 fDaughter[0].SetGood();
419 fDaughter[0].SetRefMC(daughter);
420 fDaughter[0].SetLabel(index[i]);
421 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
422 fDaughter[1].SetGood();
423 fDaughter[1].SetRefMC(daughter);
424 fDaughter[1].SetLabel(index[i]);
428 // return success if both daughters were assigned
429 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
432 fDaughter[0].Reset();
433 fDaughter[1].Reset();
438 //_____________________________________________________________________________
439 Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
442 // Loop on event and fill containers
445 // check presence of MC reference
446 if (!rsn->GetRefMC()) {
447 AliError("Need a MC to compute efficiency");
452 // must be ESD or AOD, and then use a bool to know in the rest
453 if (!rsn->IsESD() && !rsn->IsAOD()) {
454 AliError("Need to process ESD or AOD input");
458 // retrieve the MC primary vertex position
459 // and do some additional coherence checks
460 //Double_t fVertex[3] = {0.0, 0.0, 0.0};
463 //TArrayF mcVertex(3);
464 //AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
465 //genEH->PrimaryVertex(mcVertex);
466 //for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
467 npart = rsn->GetRefMCESD()->GetNumberOfTracks();
469 //for (i = 0; i < 3; i++) fVertex[i] = 0.0;
470 AliAODEvent *aod = rsn->GetRefMCAOD();
471 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
472 if (listAOD) npart = listAOD->GetEntries();
473 //AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
474 //if (mcH) mcH->GetVertex(fVertex);
477 // check number of particles
479 AliInfo("Empty event");
484 Int_t ipart, count = 0;
485 AliRsnDaughter check;
487 TObjArrayIter next(&fOutputs);
488 AliRsnListOutput *out = 0x0;
490 // loop over particles
491 for (ipart = 0; ipart < npart; ipart++) {
492 // check i-th particle
493 if (!AssignMotherAndDaughters(rsn, ipart)) continue;
494 // if assignment was successful, for first step, use MC info
495 fDaughter[0].SetRef(fDaughter[0].GetRefMC());
496 fDaughter[1].SetRef(fDaughter[1].GetRefMC());
497 fMother.SetRefEvent(rsn);
498 fMother.ComputeSum(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
499 // check rapidity range
500 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
501 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
506 while ( (out = (AliRsnListOutput *)next()) ) {
507 if (out->Fill(&fMother)) count++;
508 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));