Commit PWG2 Makefile. To compile the code of the PWG2 module one has to do: make...
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnAnalysis.cxx
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 //--------------------------------------------------------------------------------------------------------