]>
Commit | Line | Data |
---|---|---|
2f1637fb | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //------------------------------------------------------------------------- | |
17 | // Class AliRsnAnalysis | |
18 | // Reconstruction and analysis of a binary resonance | |
19 | // | |
20 | // author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it) | |
21 | //------------------------------------------------------------------------- | |
22 | ||
23 | #include <Riostream.h> | |
24 | ||
25 | #include <TH1.h> | |
26 | #include <TTree.h> | |
27 | #include <TRefArray.h> | |
28 | #include <TParticle.h> | |
29 | #include <TObjectTable.h> | |
30 | #include <TDatabasePDG.h> | |
31 | ||
32 | #include "AliRsnDaughter.h" | |
33 | #include "AliRsnDaughterCut.h" | |
34 | #include "AliRsnEvent.h" | |
35 | #include "AliRsnAnalysis.h" | |
36 | ||
37 | ClassImp(AliRsnAnalysis) | |
38 | ||
39 | //-------------------------------------------------------------------------------------------------------- | |
40 | AliRsnAnalysis::AliRsnAnalysis() | |
41 | // | |
42 | // Constructor | |
43 | // Initializes all pointers and collections to NULL. | |
44 | // | |
45 | { | |
46 | fMixHistograms = 0; | |
47 | fHistograms = 0; | |
48 | fEventsTree = 0; | |
49 | fMixPairDefs = 0; | |
50 | fPairDefs = 0; | |
51 | fPairCuts = 0; | |
52 | ||
53 | Int_t i; | |
54 | for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0; | |
55 | } | |
56 | //-------------------------------------------------------------------------------------------------------- | |
57 | void AliRsnAnalysis::AddCutPair(AliRsnDaughterCut *cut) | |
58 | // | |
59 | // Add a cut on pairs. | |
60 | // This cut is global for all pairs. | |
61 | // | |
62 | { | |
63 | if (!cut->IsPairCut()) { | |
64 | Warning("AddCutPair", "This is a single cut, cannot be added"); | |
65 | return; | |
66 | } | |
67 | ||
68 | if (!fPairCuts) fPairCuts = new TObjArray(0); | |
69 | ||
70 | fPairCuts->AddLast(cut); | |
71 | } | |
72 | //-------------------------------------------------------------------------------------------------------- | |
73 | void AliRsnAnalysis::AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut) | |
74 | // | |
75 | // Add a cut on single particles. | |
76 | // This cut must be specified for each particle type. | |
77 | // | |
78 | { | |
79 | if (cut->IsPairCut()) { | |
80 | Warning("AddCutSingle", "This is a pair cut, cannot be added"); | |
81 | return; | |
82 | } | |
83 | ||
84 | if (type >= AliPID::kElectron && type <= AliPID::kProton) { | |
85 | if (!fCuts[type]) fCuts[type] = new TObjArray(0); | |
86 | fCuts[type]->AddLast(cut); | |
87 | } | |
88 | } | |
89 | //-------------------------------------------------------------------------------------------------------- | |
90 | void AliRsnAnalysis::AddMixPairDef | |
91 | (AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2) | |
92 | // | |
93 | // Adds a new pair definition to create a new histogram, | |
94 | // for the event mixing step. | |
95 | // If the pair defs array is NULL, it is initialized here. | |
96 | // | |
97 | { | |
98 | if (!fMixPairDefs) fMixPairDefs = new TObjArray(0); | |
99 | fMixPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, 0, kFALSE) ); | |
100 | } | |
101 | //-------------------------------------------------------------------------------------------------------- | |
102 | void AliRsnAnalysis::AddPairDef | |
103 | (AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Bool_t onlyTrue) | |
104 | // | |
105 | // Adds a new pair definition to create a new histogram, | |
106 | // for the signal evaluation (same event) step. | |
107 | // If the pair defs array is NULL, it is initialized here. | |
108 | // | |
109 | { | |
110 | if (!fPairDefs) fPairDefs = new TObjArray(0); | |
111 | fPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, fTrueMotherPDG, onlyTrue) ); | |
112 | } | |
113 | //-------------------------------------------------------------------------------------------------------- | |
114 | void AliRsnAnalysis::Clear(Option_t* /* option */) | |
115 | // | |
116 | // Clear heap | |
117 | // | |
118 | { | |
119 | fHistograms->Clear("C"); | |
120 | fPairDefs->Clear("C"); | |
121 | ||
122 | fMixHistograms->Clear("C"); | |
123 | fMixPairDefs->Clear("C"); | |
124 | ||
125 | Int_t i; | |
126 | for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i]->Clear("C"); | |
127 | fPairCuts->Clear("C"); | |
128 | } | |
129 | //-------------------------------------------------------------------------------------------------------- | |
130 | Stat_t AliRsnAnalysis::EventMix(Int_t nmix, Int_t multDiffMax, Double_t vzDiffMax, Bool_t compareTotMult) | |
131 | // | |
132 | // Performs event mixing. | |
133 | // It takes the array of fMixPairDefs and stores results in fMixHistograms. | |
134 | // First parameter defines how many events must be mixed with each one. | |
135 | // Events to be mixed are chosen using the other parameters: | |
136 | // | |
137 | // - multDiffMax defines the maximum allowed difference in particle multiplicity, | |
138 | // a) when last argument is kFALSE (default) the multiplicity comparison is made with respect | |
139 | // to the particles of 'type 2' in the pair | |
140 | // b) when last argument is kTRUE, the multiplicity comparison is made with respect of total | |
141 | // particle multiplicity in the events | |
142 | // | |
143 | // - vzDiffMax defines maximum allowed difference in Z coordinate of prim. vertex. | |
144 | // | |
145 | // If one wants to exchange the particle types, one has to add both combinations of particles | |
146 | // as different pair defs. | |
147 | // | |
148 | // EXAMPLE: | |
149 | // analysis->AddMixPairDef(AliRsnEvent::kPion, '+', AliRsnEvent::kKaon, '-'); | |
150 | // analysis->AddMixPairDef(AliRsnEvent::kKaon, '-', AliRsnEvent::kPion, '+'); | |
151 | // | |
152 | { | |
153 | // allocate the histograms array | |
154 | Int_t i, npairdefs = (Int_t)fMixPairDefs->GetEntries(); | |
155 | fMixHistograms = new TObjArray(npairdefs); | |
156 | TObjArray &histos = *fMixHistograms; | |
157 | AliPairDef *pd = 0; | |
158 | for (i = 0; i < npairdefs; i++) { | |
159 | pd = (AliPairDef*)fMixPairDefs->At(i); | |
160 | histos[i] = new TH1D(Form("hmix_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax); | |
161 | } | |
162 | ||
163 | // Link events branch | |
164 | AliRsnEvent *event2 = new AliRsnEvent; | |
165 | fEventsTree->SetBranchAddress("events", &event2); | |
166 | ||
167 | // define variables to store data about particles | |
168 | Stat_t nPairs = 0; | |
169 | Int_t iev, ievmix, nEvents = (Int_t)fEventsTree->GetEntries(); | |
170 | ||
171 | // loop on events | |
172 | Int_t nmixed, mult1, mult2, diffMult; | |
173 | Double_t vz1, vz2, diffVz; | |
174 | for (iev = 0; iev < nEvents; iev++) { | |
175 | ||
176 | // get event | |
177 | event2->Clear("DELETE"); | |
178 | fEventsTree->GetEntry(iev); | |
179 | ||
180 | // copy this event into 'event 1' | |
181 | AliRsnEvent *event1 = new AliRsnEvent(*event2); | |
182 | ||
183 | // get Z position of primary vertex | |
184 | vz1 = event1->GetPrimaryVertexZ(); | |
185 | ||
186 | // if total multiplicities must be used | |
187 | // it is computed here | |
188 | if (compareTotMult) { | |
189 | mult1 = event1->GetMultiplicity(kTRUE); | |
190 | } | |
191 | else { | |
192 | mult1 = 0; | |
193 | } | |
194 | ||
195 | // message | |
196 | if (iev % 10 == 0) cout << "\rMixing event " << iev << flush; | |
197 | ||
198 | // loop on pair definitions | |
199 | for (i = 0; i < npairdefs; i++) { | |
200 | ||
201 | // get pair definition | |
202 | pd = (AliPairDef*)fMixPairDefs->At(i); | |
203 | ||
204 | // get histogram reference | |
205 | TH1D *histogram = (TH1D*)histos[i]; | |
206 | ||
207 | // get multiplicity of particle type '2' in event '1' | |
208 | if (!mult1) { | |
209 | mult1 = (Int_t)event1->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries(); | |
210 | } | |
211 | ||
212 | // loop on events for mixing | |
213 | nmixed = 0; | |
214 | ievmix = iev; | |
215 | while (nmixed < nmix) { | |
216 | ||
217 | // get other event (it starts from iev + 1, and | |
218 | // loops again to 0 if reachs end of tree) | |
219 | ievmix++; | |
220 | if (ievmix >= nEvents) ievmix -= nEvents; | |
221 | ||
222 | // if event-2-mix is equal to 'iev', then | |
223 | // a complete loop has been done, and no more | |
224 | // mixing can be done, even if they are too few | |
225 | // then the loop breaks anyway | |
226 | if (ievmix == iev) break; | |
227 | ||
228 | // get other event | |
229 | event2->Clear("DELETE"); | |
230 | fEventsTree->GetEntry(ievmix); | |
231 | ||
232 | // get event stats related to event 2 | |
233 | vz2 = event2->GetPrimaryVertexZ(); | |
234 | if (compareTotMult) { | |
235 | mult2 = event2->GetMultiplicity(kTRUE); | |
236 | } | |
237 | else { | |
238 | mult2 = event2->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries(); | |
239 | } | |
240 | ||
241 | // fill histogram if checks are passed | |
242 | diffMult = TMath::Abs(mult1 - mult2); | |
243 | diffVz = TMath::Abs(vz1 - vz2); | |
244 | if (diffVz <= vzDiffMax && diffMult <= multDiffMax) { | |
245 | //cout << ievmix << " " << flush; | |
246 | nPairs += Compute(pd, histogram, event1, event2); | |
247 | nmixed++; | |
248 | } | |
249 | } | |
250 | } | |
251 | ||
252 | event1->Clear("DELETE"); | |
253 | delete event1; | |
254 | event1 = 0; | |
255 | ||
256 | } // end loop on events | |
257 | cout << endl; | |
258 | ||
259 | return nPairs; | |
260 | } | |
261 | //-------------------------------------------------------------------------------------------------------- | |
262 | Stat_t AliRsnAnalysis::Process() | |
263 | // | |
264 | // Reads the list 'fPairDefs', and builds an inv-mass histogram for each definition. | |
265 | // For each event, particle of 'type 1' are combined with particles of 'type 2' as | |
266 | // defined in each pair definition specified in the list, taking the list of both types | |
267 | // from the same event. | |
268 | // This method is supposed to evaluate the 'signal' of the resonance, containing the peak. | |
269 | // It can also be used when one wants to evaluate the 'background' with the 'wrong sign' | |
270 | // particle pairs. | |
271 | // | |
272 | { | |
273 | // allocate the histograms array in order to contain | |
274 | // as many objects as the number of pair definitionss | |
275 | Int_t i, npairdefs = (Int_t)fPairDefs->GetEntries(); | |
276 | fHistograms = new TObjArray(npairdefs); | |
277 | TObjArray &histos = *fHistograms; | |
278 | ||
279 | // loop over pair defs in order to retrieve the particle species chosen | |
280 | // which are used to set an unique name for the output histogram | |
281 | // There will be a direct correspondence between the inder of pairdef in its array | |
282 | // and the corresponding histogram in the 'fHistograms' list. | |
283 | AliPairDef *pd = 0; | |
284 | for (i = 0; i < npairdefs; i++) { | |
285 | pd = (AliPairDef*)fPairDefs->At(i); | |
286 | histos[i] = new TH1D(Form("h_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax); | |
287 | } | |
288 | ||
289 | // Link events branch | |
290 | AliRsnEvent *event = new AliRsnEvent; | |
291 | fEventsTree->SetBranchAddress("events", &event); | |
292 | ||
293 | // define variables to store data about particles | |
294 | Stat_t nPairs = 0; | |
295 | Int_t iev, nEvents = (Int_t)fEventsTree->GetEntries(); | |
296 | ||
297 | // loop on events | |
298 | for (iev = 0; iev < nEvents; iev++) { | |
299 | ||
300 | // get event | |
301 | event->Clear("DELETE"); | |
302 | fEventsTree->GetEntry(iev); | |
303 | if (iev % 1000 == 0) cout << "\rProcessing event " << iev << flush; | |
304 | ||
305 | // loop over the collection of pair defs | |
306 | for (i = 0; i < npairdefs; i++) { | |
307 | pd = (AliPairDef*)fPairDefs->At(i); | |
308 | TH1D *histogram = (TH1D*)histos[i]; | |
309 | // invoke AliPairDef method to fill histogram | |
310 | nPairs += Compute(pd, histogram, event, event); | |
311 | } | |
312 | ||
313 | } // end loop on events | |
314 | cout << endl; | |
315 | ||
316 | return nPairs; | |
317 | } | |
318 | //-------------------------------------------------------------------------------------------------------- | |
319 | void AliRsnAnalysis::WriteHistograms() | |
320 | // | |
321 | // Writes histograms in current directory | |
322 | // | |
323 | { | |
324 | TH1D *histogram; | |
325 | TObjArrayIter iter(fHistograms); | |
326 | while ( (histogram = (TH1D*)iter.Next()) ) { | |
327 | histogram->Write(); | |
328 | } | |
329 | ||
330 | if (fMixHistograms) { | |
331 | TObjArrayIter iterMix(fMixHistograms); | |
332 | while ( (histogram = (TH1D*)iterMix.Next()) ) { | |
333 | histogram->Write(); | |
334 | } | |
335 | } | |
336 | } | |
337 | //-------------------------------------------------------------------------------------------------------- | |
338 | Stat_t AliRsnAnalysis::Compute | |
339 | (AliPairDef *pd, TH1D* &histogram, AliRsnEvent *event1, AliRsnEvent *event2) | |
340 | // | |
341 | // Adds to the specified histogram the invariant mass spectrum calculated taking | |
342 | // particles of type 1 from event 1, and particles of type 2 from event 2. | |
343 | // Events can be equal (signal) or different (background with event mixing). | |
344 | // | |
345 | { | |
346 | // define two 'cursor' objects | |
347 | AliRsnDaughter *track1 = 0, *track2 = 0; | |
348 | ||
349 | // define iterators for the two collections | |
350 | if (!event1 || !event1->GetTracks(pd->GetSign1(), pd->GetParticle1())) { | |
351 | Error("Compute", "Null pointer for particle 1 collection"); | |
352 | return 0.0; | |
353 | } | |
354 | if (!event2 || !event2->GetTracks(pd->GetSign2(), pd->GetParticle2())) { | |
355 | Error("Compute", "Null pointer for particle 2 collection"); | |
356 | return 0.0; | |
357 | } | |
358 | TObjArrayIter iter1(event1->GetTracks(pd->GetSign1(), pd->GetParticle1())); | |
359 | TObjArrayIter iter2(event2->GetTracks(pd->GetSign2(), pd->GetParticle2())); | |
360 | ||
361 | // define temporary variables for better code readability | |
362 | Stat_t nPairs = 0; | |
363 | ||
364 | // loop on particle of type 1 (in event 1) | |
365 | while ( (track1 = (AliRsnDaughter*)iter1.Next()) ) { | |
366 | ||
367 | if (fRejectFakes && (track1->GetLabel() < 0)) continue; | |
368 | if (!SingleCutCheck(pd->GetParticle1(), track1)) continue; | |
369 | ||
370 | iter2.Reset(); | |
371 | ||
372 | // loop on particles of type 2 (in event 2) | |
373 | while ( (track2 = (AliRsnDaughter*)iter2.Next()) ) { | |
374 | ||
375 | if (fRejectFakes && (track2->GetLabel() < 0)) continue; | |
376 | if (event1 == event2 && track1->GetIndex() == track2->GetIndex()) continue; | |
377 | if (!SingleCutCheck(pd->GetParticle2(), track2)) continue; | |
378 | if (!PairCutCheck(track1, track2)) continue; | |
379 | ||
380 | /* | |
381 | // check | |
382 | Char_t sign1 = (track1->GetSign() > 0) ? '+' : '-'; | |
383 | Char_t sign2 = (track2->GetSign() > 0) ? '+' : '-'; | |
384 | Info("Compute", "Particle 1: PDG = %d, sign = %c --- Particle 2: PDG = %d, sign = %c", track1->GetTruePDG(), sign1, track2->GetTruePDG(), sign2); | |
385 | */ | |
386 | ||
387 | // total 4-momentum | |
388 | track1->SetMass(pd->GetMass1()); | |
389 | track2->SetMass(pd->GetMass2()); | |
390 | AliRsnDaughter sum = (*track1) + (*track2); | |
391 | ||
392 | // if the choice to get ONLY true pairs is selected, a check is made on the mothers | |
393 | if (pd->GetOnlyTrue()) { | |
394 | // the 'sum' AliRsnDaughter object is initialized with | |
395 | // the PDG code of the common mother (if it is the case) | |
396 | if (TMath::Abs(sum.GetMotherPDG()) == TMath::Abs(fTrueMotherPDG)) { | |
397 | histogram->Fill(sum.GetMass()); | |
398 | nPairs++; | |
399 | } | |
400 | } | |
401 | else { | |
402 | histogram->Fill(sum.GetMass()); | |
403 | nPairs++; | |
404 | } | |
405 | ||
406 | } // end loop 2 | |
407 | ||
408 | } // end loop 1 | |
409 | ||
410 | return nPairs; | |
411 | } | |
412 | //-------------------------------------------------------------------------------------------------------- | |
413 | Bool_t AliRsnAnalysis::SingleCutCheck(Int_t itype, AliRsnDaughter *track) | |
414 | // | |
415 | // Checks a track against single particle cuts (if defined) | |
416 | // | |
417 | { | |
418 | if (!fCuts[itype]) return kTRUE; | |
419 | ||
420 | TObjArrayIter iter(fCuts[itype]); | |
421 | AliRsnDaughterCut *cut = 0; | |
422 | while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) { | |
423 | if (!cut->Pass(track)) return kFALSE; | |
424 | } | |
425 | ||
426 | return kTRUE; | |
427 | } | |
428 | //-------------------------------------------------------------------------------------------------------- | |
429 | Bool_t AliRsnAnalysis::PairCutCheck(AliRsnDaughter *track1, AliRsnDaughter *track2) | |
430 | // | |
431 | // Checks a pair against pair cuts (if defined) | |
432 | // | |
433 | { | |
434 | if (!fPairCuts) return kTRUE; | |
435 | ||
436 | TObjArrayIter iter(fPairCuts); | |
437 | AliRsnDaughterCut *cut = 0; | |
438 | while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) { | |
439 | if (!cut->Pass(track1, track2)) return kFALSE; | |
440 | } | |
441 | ||
442 | return kTRUE; | |
443 | } | |
444 | //-------------------------------------------------------------------------------------------------------- | |
445 | AliRsnAnalysis::AliPairDef::AliPairDef | |
446 | (AliPID::EParticleType p1, Char_t sign1, | |
447 | AliPID::EParticleType p2, Char_t sign2, Int_t pdgMother, Bool_t onlyTrue) | |
448 | // | |
449 | // Constructor for nested class | |
450 | // | |
451 | { | |
452 | fOnlyTrue = onlyTrue; | |
453 | fTrueMotherPDG = 0; | |
454 | if (fOnlyTrue) fTrueMotherPDG = pdgMother; | |
455 | ||
456 | fSign1 = sign1; | |
457 | fParticle1 = p1; | |
458 | ||
459 | fSign2 = sign2; | |
460 | fParticle2 = p2; | |
461 | ||
462 | // instance a PDG database to recovery true masses of particles | |
463 | TDatabasePDG *db = TDatabasePDG::Instance(); | |
464 | Int_t pdg1 = AliPID::ParticleCode((Int_t)p1); | |
465 | Int_t pdg2 = AliPID::ParticleCode((Int_t)p2); | |
466 | fMass1 = db->GetParticle(pdg1)->Mass(); | |
467 | fMass2 = db->GetParticle(pdg2)->Mass(); | |
468 | ||
469 | // set name according to the chosen particles | |
470 | fName.Append(Form("%s(%c)_%s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2)); | |
471 | fTitle.Append(Form("Inv. mass for %s(%c) and %s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2)); | |
472 | if (onlyTrue) { | |
473 | fName.Append("_true"); | |
474 | fTitle.Append(" (true pairs)"); | |
475 | } | |
476 | } | |
477 | //-------------------------------------------------------------------------------------------------------- | |
478 | Text_t* AliRsnAnalysis::AliPairDef::ParticleName(AliPID::EParticleType part) | |
479 | // | |
480 | // [PRIVATE] | |
481 | // Returns the name of the particle in text format | |
482 | // | |
483 | { | |
484 | if (part == AliPID::kElectron) return ("E"); | |
485 | else if (part == AliPID::kMuon) return ("Mu"); | |
486 | else if (part == AliPID::kPion) return ("Pi"); | |
487 | else if (part == AliPID::kKaon) return ("K"); | |
488 | else if (part == AliPID::kProton) return ("P"); | |
489 | else { | |
490 | Warning("ParticleName", "Unrecognized value of EParticle argument"); | |
491 | return ("???"); | |
492 | } | |
493 | } | |
494 | //-------------------------------------------------------------------------------------------------------- |