Adding AliRsnAnalysis in the RESONANCES directory
[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
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
37ClassImp(AliRsnAnalysis)
38
39//--------------------------------------------------------------------------------------------------------
40AliRsnAnalysis::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//--------------------------------------------------------------------------------------------------------
57void 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//--------------------------------------------------------------------------------------------------------
73void 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//--------------------------------------------------------------------------------------------------------
90void 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//--------------------------------------------------------------------------------------------------------
102void 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//--------------------------------------------------------------------------------------------------------
114void 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//--------------------------------------------------------------------------------------------------------
130Stat_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//--------------------------------------------------------------------------------------------------------
262Stat_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//--------------------------------------------------------------------------------------------------------
319void 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//--------------------------------------------------------------------------------------------------------
338Stat_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//--------------------------------------------------------------------------------------------------------
413Bool_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//--------------------------------------------------------------------------------------------------------
429Bool_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//--------------------------------------------------------------------------------------------------------
445AliRsnAnalysis::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//--------------------------------------------------------------------------------------------------------
478Text_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//--------------------------------------------------------------------------------------------------------