Changed scripts for new TrainSetup
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnLoopEffPair.cxx
CommitLineData
c865cb1d 1//
2// Class AliRsnLoopEffPair
3//
4// Inherits from basic AliRsnLoopEff for efficiency,
5// and computed efficiencies for single-tracks
6//
7// author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
8//
9
f34f960b 10#include <Riostream.h>
11
c865cb1d 12#include "AliStack.h"
f34f960b 13#include "AliGenEventHeader.h"
14#include "AliAODMCHeader.h"
c865cb1d 15
e187bd70 16#include "AliRsnEvent.h"
c865cb1d 17#include "AliRsnCutManager.h"
18#include "AliRsnDaughterDef.h"
19#include "AliRsnPairDef.h"
20#include "AliRsnLoopEffPair.h"
21
22ClassImp(AliRsnLoopEffPair)
23
24//_____________________________________________________________________________
25AliRsnLoopEffPair::AliRsnLoopEffPair(const char *name, AliRsnPairDef *def) :
f34f960b 26 AliRsnLoopEff(name, 1),
c865cb1d 27 fDef(def),
28 fMother()
29{
30//
31// Default constructor.
32// Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
33//
34}
35
36//_____________________________________________________________________________
61f275d1 37AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair &copy) :
c865cb1d 38 AliRsnLoopEff(copy),
39 fDef(copy.fDef),
40 fMother(copy.fMother)
41{
42//
43// Copy constrtuctor.
44//
45}
46
47//_____________________________________________________________________________
61f275d1 48AliRsnLoopEffPair &AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair &copy)
c865cb1d 49{
50//
51// Assignment operator.
52// Owned data members are meaningless for this operator.
53//
54
55 AliRsnLoopEff::operator=(copy);
e6f3a909 56 if (this == &copy)
61f275d1 57 return *this;
c865cb1d 58 fDef = copy.fDef;
59
60 return (*this);
61}
62
f34f960b 63//__________________________________________________________________________________________________
64Bool_t AliRsnLoopEffPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
65{
66//
67// Calls the appropriate assignment method
68//
69
70 // setup pointers
71 fMother.SetDaughter(0, &fDaughter[0]);
72 fMother.SetDaughter(1, &fDaughter[1]);
73 fMother.SetRefEvent(rsnEvent);
74 fDaughter[0].SetOwnerEvent(rsnEvent);
75 fDaughter[1].SetOwnerEvent(rsnEvent);
76
77 if (rsnEvent->IsESD())
78 return AssignMotherAndDaughtersESD(rsnEvent, ipart);
79 else if (rsnEvent->IsAOD())
80 return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
81 else {
82 AliError("Unrecognized input event");
83 return kFALSE;
84 }
85}
86
87//__________________________________________________________________________________________________
88Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
89{
90//
91// Gets a particle in the MC event and try to assign it to the mother.
92// If it has two daughters, retrieve them and assign also them.
93// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
94// for ESD and AOD, if available
95// ---
96// Implementation for ESD inputs
97//
98
99 AliMCEvent *mc = rsnEvent->GetRefMCESD();
100 AliStack *stack = mc->Stack();
61f275d1 101 AliMCParticle *mother = (AliMCParticle *)mc->GetTrack(ipart);
f34f960b 102 TParticle *motherP = mother->Particle();
103 Int_t ntracks = stack->GetNtrack();
61f275d1 104
f34f960b 105 // check PDG code and exit if it is wrong
106 if (TMath::Abs(motherP->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
61f275d1 107
f34f960b 108 // check number of daughters and exit if it is not 2
109 if (motherP->GetNDaughters() < 2) return kFALSE;
61f275d1 110
f34f960b 111 // check distance from primary vertex
112 TLorentzVector vprod;
113 motherP->ProductionVertex(vprod);
114 if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
115 AliDebugClass(1, "Distant production vertex");
116 return kFALSE;
117 }
61f275d1 118
f34f960b 119 // get the daughters and check their PDG code and charge:
61f275d1 120 // if they match one of the pair daughter definitions,
f34f960b 121 // assign them as MC reference of the 'fDaughter' objects
122 fDaughter[0].Reset();
123 fDaughter[1].Reset();
124 Int_t index[2] = {motherP->GetFirstDaughter(), motherP->GetLastDaughter()};
125 Int_t i, pdg;
126 Short_t charge;
127 AliMCParticle *daughter = 0x0;
128 for (i = 0; i < 2; i++) {
129 // check index for stack
130 if (index[i] < 0 || index[i] > ntracks) {
131 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
132 return kFALSE;
133 }
134 // get daughter and its PDG and charge
61f275d1 135 daughter = (AliMCParticle *)mc->GetTrack(index[i]);
f34f960b 136 pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
137 charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
138 // check if it matches one definition
139 if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
140 fDaughter[0].SetGood();
141 fDaughter[0].SetRefMC(daughter);
142 fDaughter[0].SetLabel(index[i]);
143 } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
144 fDaughter[1].SetGood();
145 fDaughter[1].SetRefMC(daughter);
146 fDaughter[1].SetLabel(index[i]);
147 }
148 }
61f275d1 149
f34f960b 150 // return success if both daughters were assigned
151 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
152 return kTRUE;
153 } else {
154 fDaughter[0].Reset();
155 fDaughter[1].Reset();
156 return kFALSE;
157 }
158}
159
160//__________________________________________________________________________________________________
161Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
162{
b63357a0 163//
164// Gets a particle in the MC event and try to assign it to the mother.
165// If it has two daughters, retrieve them and assign also them.
166// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
167// for ESD and AOD, if available
168// ---
169// Implementation for AOD inputs
170//
171
f34f960b 172 AliAODEvent *aod = rsnEvent->GetRefAOD();
61f275d1 173 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
174 AliAODMCParticle *mother = (AliAODMCParticle *)listAOD->At(ipart);
f34f960b 175 Int_t ntracks = listAOD->GetEntries();
61f275d1 176
f34f960b 177 // check PDG code and exit if it is wrong
178 if (TMath::Abs(mother->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
61f275d1 179
f34f960b 180 // check number of daughters and exit if it is not 2
181 if (mother->GetNDaughters() < 2) return kFALSE;
61f275d1 182
f34f960b 183 // check distance from primary vertex
184 Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
185 Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
186 if (dv > fMaxDistPV) {
187 AliDebugClass(1, "Distant production vertex");
188 return kFALSE;
189 }
61f275d1 190
f34f960b 191 // get the daughters and check their PDG code and charge:
61f275d1 192 // if they match one of the pair daughter definitions,
f34f960b 193 // assign them as MC reference of the 'fDaughter' objects
194 fDaughter[0].Reset();
195 fDaughter[1].Reset();
196 Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
197 Int_t i, pdg;
198 Short_t charge;
199 AliAODMCParticle *daughter = 0x0;
200 for (i = 0; i < 2; i++) {
201 // check index for stack
202 if (index[i] < 0 || index[i] > ntracks) {
203 AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
204 return kFALSE;
205 }
206 // get daughter and its PDG and charge
61f275d1 207 daughter = (AliAODMCParticle *)listAOD->At(index[i]);
f34f960b 208 pdg = TMath::Abs(daughter->GetPdgCode());
209 charge = (Short_t)(daughter->Charge() / 3);
210 // check if it matches one definition
211 if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
212 fDaughter[0].SetGood();
213 fDaughter[0].SetRefMC(daughter);
214 fDaughter[0].SetLabel(index[i]);
215 } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
216 fDaughter[1].SetGood();
217 fDaughter[1].SetRefMC(daughter);
218 fDaughter[1].SetLabel(index[i]);
219 }
220 }
61f275d1 221
f34f960b 222 // return success if both daughters were assigned
223 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
224 return kTRUE;
225 } else {
226 fDaughter[0].Reset();
227 fDaughter[1].Reset();
228 return kFALSE;
229 }
230}
231
c865cb1d 232//_____________________________________________________________________________
61f275d1 233Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent *rsn, AliRsnDaughterSelector *, AliRsnEvent *, AliRsnDaughterSelector *)
c865cb1d 234{
235//
236// Loop on event and fill containers
237//
238
239 // check event cuts
240 if (!OkEvent(rsn)) return 0;
61f275d1 241
f34f960b 242 // retrieve output
61f275d1 243 fOutput = (AliRsnListOutput *)fOutputs[0];
244
c865cb1d 245 // check presence of MC reference
246 if (!rsn->GetRefMC()) {
247 AliError("Need a MC to compute efficiency");
248 return 0;
249 }
61f275d1 250
f34f960b 251 // check presence of event
252 if (!rsn->GetRef()) {
253 AliError("Need an event to compute efficiency");
254 return 0;
255 }
61f275d1 256
c865cb1d 257 // check event type:
258 // must be ESD or AOD, and then use a bool to know in the rest
259 if (!rsn->IsESD() && !rsn->IsAOD()) {
260 AliError("Need to process ESD or AOD input");
261 return 0;
262 }
61f275d1 263
f34f960b 264 // retrieve the MC primary vertex position
265 // and do some additional coherence checks
266 Int_t i, npart = 0;
e187bd70 267 if (rsn->IsESD()) {
f34f960b 268 TArrayF mcVertex(3);
269 AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
270 genEH->PrimaryVertex(mcVertex);
271 for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
272 npart = rsn->GetRefMCESD()->GetNumberOfTracks();
c865cb1d 273 } else {
f34f960b 274 for (i = 0; i < 3; i++) fVertex[i] = 0.0;
275 AliAODEvent *aod = rsn->GetRefMCAOD();
61f275d1 276 TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
f34f960b 277 if (listAOD) npart = listAOD->GetEntries();
61f275d1 278 AliAODMCHeader *mcH = static_cast<AliAODMCHeader *>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
f34f960b 279 if (mcH) mcH->GetVertex(fVertex);
c865cb1d 280 }
61f275d1 281
c865cb1d 282 // check number of particles
283 if (!npart) {
284 AliInfo("Empty event");
285 return 0;
286 }
61f275d1 287
c865cb1d 288 // utility variables
f34f960b 289 Int_t ipart, istep, count = 0, nsteps = fSteps.GetEntries();
290 Int_t ntracks = rsn->GetAbsoluteSum();
291 AliRsnDaughter check;
61f275d1 292
c865cb1d 293 // loop over particles
294 for (ipart = 0; ipart < npart; ipart++) {
f34f960b 295 // check i-th particle
296 if (!AssignMotherAndDaughters(rsn, ipart)) continue;
297 // if assignment was successful, for first step, use MC info
298 fDaughter[0].SetRef(fDaughter[0].GetRefMC());
299 fDaughter[1].SetRef(fDaughter[1].GetRefMC());
300 fMother.SetRefEvent(rsn);
301 fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
302 fOutput->Fill(&fMother, 0);
c865cb1d 303 count++;
f34f960b 304 // for each further step, try to find two tracks which pass the related cuts
305 for (istep = 0; istep < nsteps; istep++) {
61f275d1 306 AliRsnCutManager *cuts = (AliRsnCutManager *)fSteps[istep];
f34f960b 307 fDaughter[0].SetBad();
308 fDaughter[1].SetBad();
309 for (i = 0; i < ntracks; i++) {
310 rsn->SetDaughter(check, i);
311 if (!cuts->PassCommonDaughterCuts(&check)) continue;
312 if (TMath::Abs(check.GetLabel()) == fDaughter[0].GetLabel()) {
313 if (!cuts->PassDaughter1Cuts(&check)) continue;
314 fDaughter[0].SetRef(check.GetRef());
315 fDaughter[0].SetGood();
316 } else if (TMath::Abs(check.GetLabel()) == fDaughter[1].GetLabel()) {
317 if (!cuts->PassDaughter2Cuts(&check)) continue;
318 fDaughter[1].SetRef(check.GetRef());
319 fDaughter[1].SetGood();
320 }
321 if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
322 fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
323 if (cuts->PassMotherCuts(&fMother)) {
324 fOutput->Fill(&fMother, 1 + istep);
325 break;
326 }
327 }
c865cb1d 328 }
329 }
330 }
61f275d1 331
c865cb1d 332 return count;
333}