Coverity fix
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnLoopPair.cxx
CommitLineData
c865cb1d 1//
2// *** Class AliRsnLoopPair ***
3//
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)
10//
11// authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
12// M. Vala (email: martin.vala@cern.ch)
13//
14
f34f960b 15#include <Riostream.h>
c865cb1d 16#include <TList.h>
17#include <TEntryList.h>
18
19#include "AliLog.h"
20
b63357a0 21#include "AliRsnEvent.h"
22#include "AliRsnPairDef.h"
c865cb1d 23#include "AliRsnMother.h"
24#include "AliRsnCutSet.h"
25#include "AliRsnDaughterSelector.h"
26
27#include "AliRsnLoopPair.h"
28
29ClassImp(AliRsnLoopPair)
30
31//_____________________________________________________________________________
32AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMixed) :
33 AliRsnLoop(name, isMixed),
b63357a0 34 fTrueMC(kFALSE),
c865cb1d 35 fOnlyTrue(kFALSE),
36 fCheckDecay(kFALSE),
b63357a0 37 fRangeY(1E20),
c865cb1d 38 fPairDef(def),
39 fPairCuts(0x0),
40 fMother()
41{
42//
43// Default constructor
44//
45
46 fListID[0] = -1;
47 fListID[1] = -1;
48}
49
50//_____________________________________________________________________________
51AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair& copy) :
52 AliRsnLoop(copy),
b63357a0 53 fTrueMC(copy.fTrueMC),
c865cb1d 54 fOnlyTrue(copy.fOnlyTrue),
55 fCheckDecay(copy.fCheckDecay),
b63357a0 56 fRangeY(copy.fRangeY),
c865cb1d 57 fPairDef(copy.fPairDef),
58 fPairCuts(copy.fPairCuts),
59 fMother(copy.fMother)
60{
61//
62// Copy constructor
63//
64
65 fListID[0] = copy.fListID[0];
66 fListID[1] = copy.fListID[1];
67}
68
69//_____________________________________________________________________________
70AliRsnLoopPair& AliRsnLoopPair::operator=(const AliRsnLoopPair& copy)
71{
b63357a0 72//
73// Assignment operator
74//
75
c865cb1d 76 AliRsnLoop::operator=(copy);
77
b63357a0 78 fTrueMC = copy.fTrueMC;
c865cb1d 79 fOnlyTrue = copy.fOnlyTrue;
80 fCheckDecay = copy.fCheckDecay;
b63357a0 81 fRangeY = copy.fRangeY;
c865cb1d 82 fPairDef = copy.fPairDef;
83 fPairCuts = copy.fPairCuts;
84 fMother = copy.fMother;
85 fListID[0] = copy.fListID[0];
86 fListID[1] = copy.fListID[1];
87
88 return (*this);
89}
90
91//_____________________________________________________________________________
92AliRsnLoopPair::~AliRsnLoopPair()
93{
94//
95// Destructor
96//
97}
98
99//_____________________________________________________________________________
100void AliRsnLoopPair::Print(Option_t* /*option*/) const
101{
102//
103// Prints info about pair
104//
105
106 AliRsnLoop::Print();
107}
108
109//_____________________________________________________________________________
110Bool_t AliRsnLoopPair::Init(const char *prefix, TList *list)
111{
112//
113// Initialization function.
114// Loops on all functions and eventual the ntuple, to initialize output objects.
115//
116
117 // assign data members relations
118 fMother.SetDaughter(0, &fDaughter[0]);
119 fMother.SetDaughter(1, &fDaughter[1]);
120 AliInfo(Form("[%s] Initialization", GetName()));
121
122 TString name(prefix);
123 name += '_';
124 name += GetName();
125 if (IsMixed()) name.Prepend("mix_");
126
127 return AliRsnLoop::Init(name.Data(), list);
128}
129
130//_____________________________________________________________________________
131Int_t AliRsnLoopPair::DoLoop
132(AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
133{
134//
135// Loop function.
136// Computes what is needed from passed events.
137// Returns the number of pairs successfully processed.
138//
139
140 if (fIsMixed) {
b63357a0 141 AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
c865cb1d 142 if (!evMix || !selMix) {
143 AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
144 return 0;
145 }
146 } else {
b63357a0 147 AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
c865cb1d 148 evMix = evMain;
149 selMix = selMain;
150 }
151 fMother.SetRefEvent(evMain);
152
153 // check cuts
154 if (!OkEvent(evMain)) {
b63357a0 155 AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
c865cb1d 156 return 0;
157 }
158 if (!OkEvent(evMix)) {
b63357a0 159 AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
c865cb1d 160 return 0;
161 }
162
b63357a0 163 // if it is required to loop over True MC, do this here and skip the rest of the method
164 if (fTrueMC) return LoopTrueMC(evMain);
165
c865cb1d 166 Int_t i0, i1, start, npairs = 0;
167
168 TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
169 TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
170 if (!list0 || !list1) {
171 AliError("Can't process NULL lists");
172 return 0;
173 }
b63357a0 174 AliDebugClass(3, Form("[%s]: list counts: %d, %d", GetName(), (Int_t)list0->GetN(), (Int_t)list1->GetN()));
c865cb1d 175 if (!list0->GetN() || !list1->GetN()) {
b63357a0 176 AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
c865cb1d 177 return 0;
178 }
179
180 TObjArrayIter next(&fOutputs);
181 AliRsnListOutput *out = 0x0;
182
183 for (i0 = 0; i0 < list0->GetN(); i0++) {
f34f960b 184 evMain->SetDaughter(fDaughter[0], (Int_t)list0->GetEntry(i0));
185 fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
c865cb1d 186 start = 0;
c865cb1d 187 if (!fIsMixed && list0 == list1) start = i0 + 1;
188 for (i1 = start; i1 < list1->GetN(); i1++) {
f34f960b 189 AliDebugClass(4, Form("Checking entries pair: %d (%d) with %d (%d)", (Int_t)i0, (Int_t)list0->GetEntry(i0), (Int_t)i1, (Int_t)list1->GetEntry(i1)));
190 evMix->SetDaughter(fDaughter[1], (Int_t)list1->GetEntry(i1));
191 fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
192 fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
193 fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
194 fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
195 fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
b63357a0 196 // check rapidity range
197 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
198 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
199 continue;
c865cb1d 200 }
201 // check mother
f34f960b 202 if (fOnlyTrue) {
203 if (!IsTrueMother()) {
204 AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
205 continue;
206 }
c865cb1d 207 }
b63357a0 208 // check cuts
209 if (fPairCuts) {
210 if (!fPairCuts->IsSelected(&fMother)) {
211 AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
212 continue;
213 }
214 }
c865cb1d 215 // fill outputs
216 next.Reset();
217 while ( (out = (AliRsnListOutput*)next()) ) {
218 if (out->Fill(&fMother)) npairs++;
219 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
220 }
221 }
222 }
f34f960b 223
c865cb1d 224 return npairs;
225}
226
227//_____________________________________________________________________________
f34f960b 228Bool_t AliRsnLoopPair::IsTrueMother()
c865cb1d 229{
230//
f34f960b 231// Checks to see if the mother comes from a true resonance.
232// It is triggered by the 'SetOnlyTrue()' function
c865cb1d 233//
f34f960b 234
235 // check #1:
236 // daughters have same mother with the right PDG code
237 Int_t commonPDG = fMother.CommonMother();
238 if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
b63357a0 239 AliDebugClass(1, "Found a true mother");
c865cb1d 240
f34f960b 241 // check #2:
242 // checks if daughter have the right particle type
243 // (activated by fCheckDecay)
244 if (fCheckDecay) {
245 AliRsnDaughterDef &def1 = fPairDef->GetDef1();
246 AliRsnDaughterDef &def2 = fPairDef->GetDef2();
247 if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
248 if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
c865cb1d 249 }
b63357a0 250 AliDebugClass(1, "Decay products match");
c865cb1d 251
f34f960b 252 return kTRUE;
c865cb1d 253}
b63357a0 254
255//__________________________________________________________________________________________________
256Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
257{
258//
259// Calls the appropriate assignment method
260//
261
262 // setup pointers
263 fMother.SetDaughter(0, &fDaughter[0]);
264 fMother.SetDaughter(1, &fDaughter[1]);
265 fMother.SetRefEvent(rsnEvent);
266 fDaughter[0].SetOwnerEvent(rsnEvent);
267 fDaughter[1].SetOwnerEvent(rsnEvent);
268
269 if (rsnEvent->IsESD())
270 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
271 else if (rsnEvent->IsAOD())
272 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
273 else {
274 AliError("Unrecognized input event");
275 return kFALSE;
276 }
277}
278
279//__________________________________________________________________________________________________
280Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
281{
282//
283// Gets a particle in the MC event and try to assign it to the mother.
284// If it has two daughters, retrieve them and assign also them.
285// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
286// for ESD and AOD, if available
287// ---
288// Implementation for ESD inputs
289//
290
291 AliMCEvent *mc = rsnEvent->GetRefMCESD();
292 AliStack *stack = mc->Stack();
293 AliMCParticle *mother = (AliMCParticle*)mc->GetTrack(ipart);
294 TParticle *motherP = mother->Particle();
295 Int_t ntracks = stack->GetNtrack();
296
297 // check PDG code and exit if it is wrong
298 if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
299
300 // check number of daughters and exit if it is not 2
301 if (motherP->GetNDaughters() < 2) return kFALSE;
302 //if (!stack->IsPhysicalPrimary(ipart)) return kFALSE;
303
304 /*
305 // check distance from primary vertex
306 TLorentzVector vprod;
307 motherP->ProductionVertex(vprod);
308 if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
309 AliDebugClass(1, "Distant production vertex");
310 return kFALSE;
311 }
312 */
313
314 // get the daughters and check their PDG code and charge:
315 // if they match one of the pair daughter definitions,
316 // assign them as MC reference of the 'fDaughter' objects
317 fDaughter[0].Reset();
318 fDaughter[1].Reset();
319 Int_t index[2] = {motherP->GetDaughter(0), motherP->GetDaughter(1)};
320 Int_t i, pdg;
321 Short_t charge;
322 AliMCParticle *daughter = 0x0;
323 for (i = 0; i < 2; i++) {
324 // check index for stack
325 if (index[i] < 0 || index[i] > ntracks) {
326 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
327 return kFALSE;
328 }
329 // get daughter and its PDG and charge
330 daughter = (AliMCParticle*)mc->GetTrack(index[i]);
331 pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
332 charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
333 // check if it matches one definition
334 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
335 fDaughter[0].SetGood();
336 fDaughter[0].SetRefMC(daughter);
337 fDaughter[0].SetLabel(index[i]);
338 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
339 fDaughter[1].SetGood();
340 fDaughter[1].SetRefMC(daughter);
341 fDaughter[1].SetLabel(index[i]);
342 }
343 }
344
345 // return success if both daughters were assigned
346 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
347 return kTRUE;
348 } else {
349 fDaughter[0].Reset();
350 fDaughter[1].Reset();
351 return kFALSE;
352 }
353}
354
355//__________________________________________________________________________________________________
356Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
357{
358//
359// Gets a particle in the MC event and try to assign it to the mother.
360// If it has two daughters, retrieve them and assign also them.
361// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
362// for ESD and AOD, if available
363// ---
364// Implementation for AOD inputs
365//
366
367 AliAODEvent *aod = rsnEvent->GetRefAOD();
368 TClonesArray *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
369 AliAODMCParticle *mother = (AliAODMCParticle*)listAOD->At(ipart);
370
371 // check PDG code and exit if it is wrong
372 if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
373
374 // check number of daughters and exit if it is not 2
375 if (mother->GetNDaughters() < 2) return kFALSE;
376 if (!mother->IsPrimary()) return kFALSE;
377
378 /*
379 // check distance from primary vertex
380 Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
381 Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
382 if (dv > fMaxDistPV) {
383 AliDebugClass(1, "Distant production vertex");
384 return kFALSE;
385 }
386 */
387
388 // get the daughters and check their PDG code and charge:
389 // if they match one of the pair daughter definitions,
390 // assign them as MC reference of the 'fDaughter' objects
391 fDaughter[0].Reset();
392 fDaughter[1].Reset();
393 Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
394 Int_t ntracks = listAOD->GetEntriesFast();
395 Int_t i, pdg;
396 Short_t charge;
397 AliAODMCParticle *daughter = 0x0;
398 for (i = 0; i < 2; i++) {
399 // check index for stack
400 if (index[i] < 0 || index[i] > ntracks) {
401 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
402 return kFALSE;
403 }
404 // get daughter and its PDG and charge
405 daughter = (AliAODMCParticle*)listAOD->At(index[i]);
406 pdg = TMath::Abs(daughter->GetPdgCode());
407 charge = (Short_t)(daughter->Charge() / 3);
408 // check if it matches one definition
409 if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
410 fDaughter[0].SetGood();
411 fDaughter[0].SetRefMC(daughter);
412 fDaughter[0].SetLabel(index[i]);
413 } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
414 fDaughter[1].SetGood();
415 fDaughter[1].SetRefMC(daughter);
416 fDaughter[1].SetLabel(index[i]);
417 }
418 }
419
420 // return success if both daughters were assigned
421 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
422 return kTRUE;
423 } else {
424 fDaughter[0].Reset();
425 fDaughter[1].Reset();
426 return kFALSE;
427 }
428}
429
430//_____________________________________________________________________________
431Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
432{
433//
434// Loop on event and fill containers
435//
436
437 // check presence of MC reference
438 if (!rsn->GetRefMC()) {
439 AliError("Need a MC to compute efficiency");
440 return 0;
441 }
442
443 // check event type:
444 // must be ESD or AOD, and then use a bool to know in the rest
445 if (!rsn->IsESD() && !rsn->IsAOD()) {
446 AliError("Need to process ESD or AOD input");
447 return 0;
448 }
449
450 // retrieve the MC primary vertex position
451 // and do some additional coherence checks
452 //Double_t fVertex[3] = {0.0, 0.0, 0.0};
453 Int_t npart = 0;
454 if (rsn->IsESD()) {
455 //TArrayF mcVertex(3);
456 //AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
457 //genEH->PrimaryVertex(mcVertex);
458 //for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
459 npart = rsn->GetRefMCESD()->GetNumberOfTracks();
460 } else {
461 //for (i = 0; i < 3; i++) fVertex[i] = 0.0;
462 AliAODEvent *aod = rsn->GetRefMCAOD();
463 TClonesArray *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
464 if (listAOD) npart = listAOD->GetEntries();
465 //AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
466 //if (mcH) mcH->GetVertex(fVertex);
467 }
468
469 // check number of particles
470 if (!npart) {
471 AliInfo("Empty event");
472 return 0;
473 }
474
475 // utility variables
476 Int_t ipart, count = 0;
477 AliRsnDaughter check;
478
479 TObjArrayIter next(&fOutputs);
480 AliRsnListOutput *out = 0x0;
481
482 // loop over particles
483 for (ipart = 0; ipart < npart; ipart++) {
484 // check i-th particle
485 if (!AssignMotherAndDaughters(rsn, ipart)) continue;
486 // if assignment was successful, for first step, use MC info
487 fDaughter[0].SetRef(fDaughter[0].GetRefMC());
488 fDaughter[1].SetRef(fDaughter[1].GetRefMC());
489 fMother.SetRefEvent(rsn);
490 fMother.ComputeSum(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
491 // check rapidity range
492 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
493 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
494 continue;
495 }
496 // fill outputs
497 next.Reset();
498 while ( (out = (AliRsnListOutput*)next()) ) {
499 if (out->Fill(&fMother)) count++;
500 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
501 }
502 }
503
504 return count;
505}