example macros to run on proof
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnLoopPair.cxx
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
15 #include <Riostream.h>
16 #include <TList.h>
17 #include <TEntryList.h>
18
19 #include "AliLog.h"
20
21 #include "AliRsnEvent.h"
22 #include "AliRsnPairDef.h"
23 #include "AliRsnMother.h"
24 #include "AliRsnCutSet.h"
25 #include "AliRsnDaughterSelector.h"
26
27 #include "AliRsnLoopPair.h"
28
29 ClassImp(AliRsnLoopPair)
30
31 //_____________________________________________________________________________
32 AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMixed) :
33    AliRsnLoop(name, isMixed),
34    fTrueMC(kFALSE),
35    fOnlyTrue(kFALSE),
36    fCheckDecay(kFALSE),
37    fRangeY(1E20),
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 //_____________________________________________________________________________
51 AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair& copy) :
52    AliRsnLoop(copy),
53    fTrueMC(copy.fTrueMC),
54    fOnlyTrue(copy.fOnlyTrue),
55    fCheckDecay(copy.fCheckDecay),
56    fRangeY(copy.fRangeY),
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 //_____________________________________________________________________________
70 AliRsnLoopPair& AliRsnLoopPair::operator=(const AliRsnLoopPair& copy)
71 {
72 //
73 // Assignment operator
74 //
75
76    AliRsnLoop::operator=(copy);
77    
78    fTrueMC = copy.fTrueMC;
79    fOnlyTrue = copy.fOnlyTrue;
80    fCheckDecay = copy.fCheckDecay;
81    fRangeY = copy.fRangeY;
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 //_____________________________________________________________________________
92 AliRsnLoopPair::~AliRsnLoopPair()
93 {
94 //
95 // Destructor
96 //
97 }
98
99 //_____________________________________________________________________________
100 void AliRsnLoopPair::Print(Option_t* /*option*/) const
101 {
102 //
103 // Prints info about pair
104 //
105
106    AliRsnLoop::Print();
107 }
108
109 //_____________________________________________________________________________
110 Bool_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 //_____________________________________________________________________________
131 Int_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) {
141       AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
142       if (!evMix || !selMix) {
143          AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
144          return 0;
145       }
146    } else {
147       AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
148       evMix = evMain;
149       selMix = selMain;
150    }
151    fMother.SetRefEvent(evMain);
152    
153    // check cuts
154    if (!OkEvent(evMain)) {
155       AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
156       return 0;
157    }
158    if (!OkEvent(evMix)) {
159       AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
160       return 0;
161    }
162    
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    
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    }
174    AliDebugClass(3, Form("[%s]: list counts: %d, %d", GetName(), (Int_t)list0->GetN(), (Int_t)list1->GetN()));
175    if (!list0->GetN() || !list1->GetN()) {
176       AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
177       return 0;
178    }
179    
180    TObjArrayIter next(&fOutputs);
181    AliRsnListOutput *out = 0x0;
182    
183    for (i0 = 0; i0 < list0->GetN(); i0++) {
184       evMain->SetDaughter(fDaughter[0], (Int_t)list0->GetEntry(i0));
185       fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
186       start = 0;
187       if (!fIsMixed && list0 == list1) start = i0 + 1;
188       for (i1 = start; i1 < list1->GetN(); i1++) {
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());
196          // check rapidity range
197          if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
198             AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
199             continue;
200          }
201          // check mother
202          if (fOnlyTrue) {
203             if (!IsTrueMother()) {
204                AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
205                continue;
206             }
207          }
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          }
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    }
223
224    return npairs;
225 }
226
227 //_____________________________________________________________________________
228 Bool_t AliRsnLoopPair::IsTrueMother()
229 {
230 //
231 // Checks to see if the mother comes from a true resonance.
232 // It is triggered by the 'SetOnlyTrue()' function
233 //
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;
239    AliDebugClass(1, "Found a true mother");
240    
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;
249    }
250    AliDebugClass(1, "Decay products match");
251    
252    return kTRUE;
253 }
254
255 //__________________________________________________________________________________________________
256 Bool_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 //__________________________________________________________________________________________________
280 Bool_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 //__________________________________________________________________________________________________
356 Bool_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 //_____________________________________________________________________________
431 Int_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 }