Changed scripts for new TrainSetup
[u/mrichter/AliRoot.git] / PWGLF / 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),
52e6652d 36 fUseMCRef(kFALSE),
c865cb1d 37 fCheckDecay(kFALSE),
b63357a0 38 fRangeY(1E20),
c865cb1d 39 fPairDef(def),
40 fPairCuts(0x0),
41 fMother()
42{
43//
44// Default constructor
45//
46
47 fListID[0] = -1;
48 fListID[1] = -1;
49}
50
51//_____________________________________________________________________________
52e6652d 52AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair &copy) :
c865cb1d 53 AliRsnLoop(copy),
b63357a0 54 fTrueMC(copy.fTrueMC),
c865cb1d 55 fOnlyTrue(copy.fOnlyTrue),
52e6652d 56 fUseMCRef(copy.fUseMCRef),
c865cb1d 57 fCheckDecay(copy.fCheckDecay),
b63357a0 58 fRangeY(copy.fRangeY),
c865cb1d 59 fPairDef(copy.fPairDef),
60 fPairCuts(copy.fPairCuts),
61 fMother(copy.fMother)
62{
63//
64// Copy constructor
65//
66
67 fListID[0] = copy.fListID[0];
68 fListID[1] = copy.fListID[1];
69}
70
71//_____________________________________________________________________________
52e6652d 72AliRsnLoopPair &AliRsnLoopPair::operator=(const AliRsnLoopPair &copy)
c865cb1d 73{
b63357a0 74//
75// Assignment operator
76//
77
c865cb1d 78 AliRsnLoop::operator=(copy);
e6f3a909 79 if (this == &copy)
61f275d1 80 return *this;
b63357a0 81 fTrueMC = copy.fTrueMC;
c865cb1d 82 fOnlyTrue = copy.fOnlyTrue;
52e6652d 83 fUseMCRef = copy.fUseMCRef;
c865cb1d 84 fCheckDecay = copy.fCheckDecay;
b63357a0 85 fRangeY = copy.fRangeY;
c865cb1d 86 fPairDef = copy.fPairDef;
87 fPairCuts = copy.fPairCuts;
88 fMother = copy.fMother;
89 fListID[0] = copy.fListID[0];
90 fListID[1] = copy.fListID[1];
91
92 return (*this);
93}
94
95//_____________________________________________________________________________
96AliRsnLoopPair::~AliRsnLoopPair()
97{
98//
99// Destructor
100//
101}
102
103//_____________________________________________________________________________
3da8cef7 104void AliRsnLoopPair::Print(Option_t * /*option*/) const
c865cb1d 105{
106//
107// Prints info about pair
108//
109
110 AliRsnLoop::Print();
111}
112
113//_____________________________________________________________________________
114Bool_t AliRsnLoopPair::Init(const char *prefix, TList *list)
115{
116//
117// Initialization function.
118// Loops on all functions and eventual the ntuple, to initialize output objects.
119//
120
121 // assign data members relations
122 fMother.SetDaughter(0, &fDaughter[0]);
123 fMother.SetDaughter(1, &fDaughter[1]);
124 AliInfo(Form("[%s] Initialization", GetName()));
52e6652d 125
c865cb1d 126 TString name(prefix);
3da8cef7 127 name += '.';
c865cb1d 128 name += GetName();
3da8cef7 129// if (IsMixed()) name.Append("_mix");
52e6652d 130
c865cb1d 131 return AliRsnLoop::Init(name.Data(), list);
132}
133
134//_____________________________________________________________________________
135Int_t AliRsnLoopPair::DoLoop
136(AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
137{
138//
139// Loop function.
140// Computes what is needed from passed events.
141// Returns the number of pairs successfully processed.
142//
143
144 if (fIsMixed) {
b63357a0 145 AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
c865cb1d 146 if (!evMix || !selMix) {
147 AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
148 return 0;
149 }
150 } else {
b63357a0 151 AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
c865cb1d 152 evMix = evMain;
153 selMix = selMain;
154 }
155 fMother.SetRefEvent(evMain);
52e6652d 156
c865cb1d 157 // check cuts
158 if (!OkEvent(evMain)) {
b63357a0 159 AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
c865cb1d 160 return 0;
161 }
162 if (!OkEvent(evMix)) {
b63357a0 163 AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
c865cb1d 164 return 0;
165 }
52e6652d 166
b63357a0 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);
52e6652d 169
c865cb1d 170 Int_t i0, i1, start, npairs = 0;
52e6652d 171
c865cb1d 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");
176 return 0;
177 }
52e6652d 178 AliDebugClass(3, Form("[%s]: list counts: %lld, %lld", GetName(), list0->GetN(), list1->GetN()));
c865cb1d 179 if (!list0->GetN() || !list1->GetN()) {
b63357a0 180 AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
c865cb1d 181 return 0;
182 }
52e6652d 183
c865cb1d 184 TObjArrayIter next(&fOutputs);
185 AliRsnListOutput *out = 0x0;
52e6652d 186 Long64_t iEntry1,iEntry2;
c865cb1d 187 for (i0 = 0; i0 < list0->GetN(); i0++) {
52e6652d 188 iEntry1 = list0->GetEntry(i0);
189 evMain->SetDaughter(fDaughter[0], iEntry1,fUseMCRef);
f34f960b 190 fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
c865cb1d 191 start = 0;
c865cb1d 192 if (!fIsMixed && list0 == list1) start = i0 + 1;
193 for (i1 = start; i1 < list1->GetN(); i1++) {
52e6652d 194 iEntry2 = list1->GetEntry(i1);
61f275d1 195 if (iEntry1 == iEntry2) continue;
52e6652d 196 AliDebugClass(4, Form("Checking entries pair: %d (%lld) with %d (%lld)", i0, iEntry1, i1, iEntry2));
197 evMix->SetDaughter(fDaughter[1], iEntry2,fUseMCRef);
f34f960b 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());
b63357a0 203 // check rapidity range
204 if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
205 AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
206 continue;
c865cb1d 207 }
208 // check mother
f34f960b 209 if (fOnlyTrue) {
210 if (!IsTrueMother()) {
211 AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
212 continue;
213 }
c865cb1d 214 }
b63357a0 215 // check cuts
216 if (fPairCuts) {
217 if (!fPairCuts->IsSelected(&fMother)) {
218 AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
219 continue;
220 }
221 }
c865cb1d 222 // fill outputs
223 next.Reset();
52e6652d 224 while ( (out = (AliRsnListOutput *)next()) ) {
c865cb1d 225 if (out->Fill(&fMother)) npairs++;
226 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
227 }
228 }
229 }
f34f960b 230
c865cb1d 231 return npairs;
232}
233
234//_____________________________________________________________________________
f34f960b 235Bool_t AliRsnLoopPair::IsTrueMother()
c865cb1d 236{
237//
f34f960b 238// Checks to see if the mother comes from a true resonance.
239// It is triggered by the 'SetOnlyTrue()' function
c865cb1d 240//
52e6652d 241
f34f960b 242 // check #1:
243 // daughters have same mother with the right PDG code
244 Int_t commonPDG = fMother.CommonMother();
245 if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
b63357a0 246 AliDebugClass(1, "Found a true mother");
52e6652d 247
f34f960b 248 // check #2:
249 // checks if daughter have the right particle type
250 // (activated by fCheckDecay)
251 if (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;
c865cb1d 256 }
b63357a0 257 AliDebugClass(1, "Decay products match");
52e6652d 258
f34f960b 259 return kTRUE;
c865cb1d 260}
b63357a0 261
262//__________________________________________________________________________________________________
263Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
264{
265//
266// Calls the appropriate assignment method
267//
268
269 // setup pointers
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);
275
276 if (rsnEvent->IsESD())
277 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
278 else if (rsnEvent->IsAOD())
279 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
280 else {
281 AliError("Unrecognized input event");
282 return kFALSE;
283 }
284}
285
286//__________________________________________________________________________________________________
287Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
288{
289//
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
294// ---
295// Implementation for ESD inputs
296//
297
298 AliMCEvent *mc = rsnEvent->GetRefMCESD();
299 AliStack *stack = mc->Stack();
52e6652d 300 AliMCParticle *mother = (AliMCParticle *)mc->GetTrack(ipart);
b63357a0 301 TParticle *motherP = mother->Particle();
302 Int_t ntracks = stack->GetNtrack();
52e6652d 303
b63357a0 304 // check PDG code and exit if it is wrong
305 if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
52e6652d 306
b63357a0 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;
52e6652d 310
b63357a0 311 /*
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");
317 return kFALSE;
318 }
319 */
52e6652d 320
b63357a0 321 // get the daughters and check their PDG code and charge:
52e6652d 322 // if they match one of the pair daughter definitions,
b63357a0 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)};
327 Int_t i, pdg;
328 Short_t charge;
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));
334 return kFALSE;
335 }
336 // get daughter and its PDG and charge
52e6652d 337 daughter = (AliMCParticle *)mc->GetTrack(index[i]);
b63357a0 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]);
349 }
350 }
52e6652d 351
b63357a0 352 // return success if both daughters were assigned
353 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
354 return kTRUE;
355 } else {
356 fDaughter[0].Reset();
357 fDaughter[1].Reset();
358 return kFALSE;
359 }
360}
361
362//__________________________________________________________________________________________________
363Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
364{
365//
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
370// ---
371// Implementation for AOD inputs
372//
373
374 AliAODEvent *aod = rsnEvent->GetRefAOD();
52e6652d 375 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
376 AliAODMCParticle *mother = (AliAODMCParticle *)listAOD->At(ipart);
377
b63357a0 378 // check PDG code and exit if it is wrong
379 if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
52e6652d 380
b63357a0 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;
52e6652d 384
b63357a0 385 /*
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");
391 return kFALSE;
392 }
393 */
52e6652d 394
b63357a0 395 // get the daughters and check their PDG code and charge:
52e6652d 396 // if they match one of the pair daughter definitions,
b63357a0 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();
402 Int_t i, pdg;
403 Short_t charge;
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));
409 return kFALSE;
410 }
411 // get daughter and its PDG and charge
52e6652d 412 daughter = (AliAODMCParticle *)listAOD->At(index[i]);
b63357a0 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]);
424 }
425 }
52e6652d 426
b63357a0 427 // return success if both daughters were assigned
428 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
429 return kTRUE;
430 } else {
431 fDaughter[0].Reset();
432 fDaughter[1].Reset();
433 return kFALSE;
434 }
435}
436
437//_____________________________________________________________________________
438Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
439{
440//
441// Loop on event and fill containers
442//
443
444 // check presence of MC reference
445 if (!rsn->GetRefMC()) {
446 AliError("Need a MC to compute efficiency");
447 return 0;
448 }
52e6652d 449
b63357a0 450 // check event type:
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");
454 return 0;
455 }
52e6652d 456
b63357a0 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};
460 Int_t npart = 0;
461 if (rsn->IsESD()) {
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();
467 } else {
468 //for (i = 0; i < 3; i++) fVertex[i] = 0.0;
469 AliAODEvent *aod = rsn->GetRefMCAOD();
52e6652d 470 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
b63357a0 471 if (listAOD) npart = listAOD->GetEntries();
472 //AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
473 //if (mcH) mcH->GetVertex(fVertex);
474 }
52e6652d 475
b63357a0 476 // check number of particles
477 if (!npart) {
478 AliInfo("Empty event");
479 return 0;
480 }
52e6652d 481
b63357a0 482 // utility variables
483 Int_t ipart, count = 0;
484 AliRsnDaughter check;
52e6652d 485
b63357a0 486 TObjArrayIter next(&fOutputs);
487 AliRsnListOutput *out = 0x0;
52e6652d 488
b63357a0 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()));
501 continue;
502 }
503 // fill outputs
504 next.Reset();
52e6652d 505 while ( (out = (AliRsnListOutput *)next()) ) {
b63357a0 506 if (out->Fill(&fMother)) count++;
507 else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
508 }
509 }
52e6652d 510
b63357a0 511 return count;
512}
52e6652d 513