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