Package upgrade.
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnReader.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 AliRsnReader ========
18 //
19 // This object reads a 'standard' event and converts it into the internal
20 // format used for resonance analysis (AliRsnEvent).
21 // 'Standard' event means ESD, standard AOD and MC event.
22 //
23 // The input-2-AliRsnEvent conversion is done through a class which reads
24 // from AliAnalysisTaskSE, which is the standard analysis object. 
25 // This class creates the AliRsnEvent's before the input event is read, 
26 // so this class has not to 'create' a new outpu event, but instead it has 
27 // to 'fill' one which has already been created elsewhere.
28 // Then, the methods provided here accept an AliRsnEvent as argument passed
29 // by reference, and they 'fill' this object using the data from the inputs
30 // passed to them.
31 // 
32 // author: A. Pulvirenti
33 // email : alberto.pulvirenti@ct.infn.it
34 //
35
36 #include <TString.h>
37
38 #include "AliLog.h"
39
40 #include "AliVEvent.h"
41
42 #include "AliESDEvent.h"
43 #include "AliESDtrack.h"
44 #include "AliESDVertex.h"
45
46 #include "AliAODEvent.h"
47 #include "AliAODTrack.h"
48 #include "AliAODVertex.h"
49
50 #include "AliStack.h"
51 #include "AliMCEvent.h"
52 #include "AliMCParticle.h"
53 #include "AliGenEventHeader.h"
54
55 #include "AliRsnMCInfo.h"
56 #include "AliRsnDaughter.h"
57 #include "AliRsnEvent.h"
58 #include "AliRsnPIDWeightsMgr.h"
59 #include "AliRsnReader.h"
60
61 ClassImp(AliRsnReader)
62
63 //_____________________________________________________________________________
64 AliRsnReader::AliRsnReader(ESource source, AliRsnPIDWeightsMgr *mgr) :
65   TObject(),
66   fSource(source),
67   fCheckSplit(kFALSE),
68   fRejectFakes(kFALSE),
69   fWeightsMgr(mgr)
70 {
71 //
72 // Constructor.
73 // Initializes the base-type data members:
74 //   - management of fake tracks
75 //
76 }
77
78 //_____________________________________________________________________________
79 AliRsnReader::AliRsnReader(const AliRsnReader &copy) : 
80   TObject(copy),
81   fSource(copy.fSource),
82   fCheckSplit(copy.fCheckSplit),
83   fRejectFakes(copy.fRejectFakes),
84   fWeightsMgr(copy.fWeightsMgr)
85 {
86 //
87 // Copy constructor.
88 //
89 }
90
91 //_____________________________________________________________________________
92 AliRsnReader& AliRsnReader::operator=(const AliRsnReader &copy)
93 {
94 //
95 // Assignment operator.
96 //
97
98     fSource = copy.fSource;
99     fCheckSplit = copy.fCheckSplit;
100     fRejectFakes = copy.fRejectFakes;
101     fWeightsMgr = copy.fWeightsMgr;
102     return (*this);
103 }
104
105 //_____________________________________________________________________________
106 Bool_t AliRsnReader::Fill(AliRsnEvent *rsn, AliVEvent *event, AliMCEvent *mc)
107 {
108 //
109 // According to the class type of event and the selected source
110 // recalls one of the private reading methods to fill the RsnEvent
111 // passed by reference as first argument.
112 //
113
114     TString str(event->ClassName());
115
116     switch (fSource) {
117         case kESD:
118             if (!str.Contains("AliESDEvent")) {
119                 AliError("A reader set to 'kESD' or 'kESDTPC' can read only ESD events");
120                 return kFALSE;
121             }
122             return FillFromESD(rsn, (AliESDEvent*)event, mc);
123         case kESDTPC:
124             if (!str.Contains("AliESDEvent")) {
125                 AliError("A reader set to 'kESD' or 'kESDTPC' can read only ESD events");
126                 return kFALSE;
127             }
128             return FillFromESD(rsn, (AliESDEvent*)event, mc, kTRUE);
129         case kAOD:
130             if (!str.Contains("AliAODEvent")) {
131                 AliError("A reader set to 'kAOD' can read only AOD events");
132                 return kFALSE;
133             }
134             return FillFromAOD(rsn, (AliAODEvent*)event, mc);
135         case kMC:
136             if (!str.Contains("AliMCEvent")) {
137                 AliError("A reader set to 'kMC' can read only MC events");
138                 return kFALSE;
139             }
140             return FillFromMC(rsn, (AliMCEvent*)event);
141         default:
142             return kFALSE;
143     }
144 }
145
146 //_____________________________________________________________________________
147 Bool_t AliRsnReader::FillFromESD
148 (AliRsnEvent *rsn, AliESDEvent *esd, AliMCEvent *mc, Bool_t useTPCOnly)
149 {
150 //
151 // Filler from an ESD event.
152 // Stores all tracks (if a filter is defined, it will store
153 // only the ones which survive the cuts).
154 // If a reference MC event is provided, it is used to store
155 // the MC informations for each track (true PDG code, 
156 // GEANT label of mother, PDG code of mother, if any).
157 // When this is used, the 'source' flag of the output
158 // AliRsnEvent object will be set to 'kESD'.
159 //
160     
161     // retrieve stack (if possible)
162     AliStack *stack = 0x0;
163     if (mc) stack = mc->Stack();
164     
165     // get number of tracks
166     Int_t ntracks = esd->GetNumberOfTracks();
167     /*
168     if (!ntracks) {
169        AliWarning("No tracks in this event");
170        return kFALSE;
171     }
172     */
173     
174     // if required with the flag, scans the event
175     // and searches all split tracks (= 2 tracks with the same label);
176     // for each pair of split tracks, only the better (best chi2) is kept
177     // and the other is rejected: this info is stored into a Boolean array
178     Int_t i1, i2, lab1, lab2;
179     Bool_t *accept = new Bool_t[ntracks];
180     for (i1 = 0; i1 < ntracks; i1++) accept[i1] = kTRUE;
181     if (fCheckSplit) {
182         for (i1 = 0; i1 < ntracks; i1++) {
183             AliESDtrack *trk1 = esd->GetTrack(i1);
184             lab1 = TMath::Abs(trk1->GetLabel());
185             for (i2 = i1+1; i2 < ntracks; i2++) {
186                 AliESDtrack *trk2 = esd->GetTrack(i2);
187                 lab2 = TMath::Abs(trk2->GetLabel());
188                 // check if labels are equal
189                 if (lab1 == lab2) {
190                     if (trk1->GetConstrainedChi2() < trk2->GetConstrainedChi2()) {
191                         accept[i1] = kTRUE;
192                         accept[i2] = kFALSE;
193                     }
194                     else {
195                         accept[i1] = kFALSE;
196                         accept[i2] = kTRUE;
197                     }
198                 }
199             }
200         }
201     }
202     
203     // get primary vertex
204     Double_t vertex[3];
205     if (!useTPCOnly) {
206         vertex[0] = esd->GetVertex()->GetXv();
207         vertex[1] = esd->GetVertex()->GetYv();
208         vertex[2] = esd->GetVertex()->GetZv();
209     }
210     else {
211         vertex[0] = esd->GetPrimaryVertexTPC()->GetXv();
212         vertex[1] = esd->GetPrimaryVertexTPC()->GetYv();
213         vertex[2] = esd->GetPrimaryVertexTPC()->GetZv();
214     }
215     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
216     
217     // store tracks from ESD
218     Int_t  i, index, label, labmum;
219     Bool_t check;
220     AliRsnDaughter temp;
221     for (index = 0; index < ntracks; index++) {
222         // skip track recognized as the worse one in a splitted pair
223         if (!accept[index]) {
224             AliInfo(Form("Rejecting split track #%d in this event", index));
225             continue;
226         }
227         // get and check ESD track
228         AliESDtrack *esdTrack = esd->GetTrack(index);
229         label = esdTrack->GetLabel();
230         if (fRejectFakes && (label < 0)) continue;
231         // copy ESD track data into RsnDaughter
232         // if unsuccessful, this track is skipped
233         check = temp.Adopt(esdTrack, useTPCOnly);
234         if (!check) continue;
235         // if the AliRsnWeightsMgr object is initialized
236         // this means that ESD PID weights are not used
237         // and they are computed according to the Weight manager
238         if (fWeightsMgr) {
239             //AliInfo("Using customized weights");
240             //AliInfo(Form("ESD weights = %f %f %f %f %f", temp.PID()[0], temp.PID()[1], temp.PID()[2], temp.PID()[3], temp.PID()[4], temp.PID()[5]));
241             esdTrack->GetITSpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kITS));
242             esdTrack->GetTPCpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTPC));
243             esdTrack->GetTRDpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTRD));
244             esdTrack->GetTOFpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kTOF));
245             esdTrack->GetHMPIDpid(fWeightsMgr->GetWeightArray(AliRsnPIDWeightsMgr::kHMPID));
246             for (i = 0; i < AliRsnPID::kSpecies; i++) {
247                 temp.SetPIDWeight(i, fWeightsMgr->GetWeight((AliRsnPID::EType)i, temp.Pt()));
248             }
249             //AliInfo(Form("Used weights = %f %f %f %f %f", temp.PID()[0], temp.PID()[1], temp.PID()[2], temp.PID()[3], temp.PID()[4], temp.PID()[5]));
250         }
251         else {
252             //AliInfo("Using standard ESD weights");
253         }
254         
255         // if stack is present, copy MC info
256         if (stack) {
257             TParticle *part = stack->Particle(TMath::Abs(label));
258             if (part) {
259                 temp.InitMCInfo(part);
260                 labmum = part->GetFirstMother();
261                 if (labmum >= 0) {
262                     TParticle *mum = stack->Particle(labmum);
263                     temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
264                 }
265             }
266         }
267         // set index and label and add this object to the output container
268         temp.SetIndex(index);
269         temp.SetLabel(label);
270         AliRsnDaughter *ptr = rsn->AddTrack(temp);
271         // if problems occurred while storing, that pointer is NULL
272         if (!ptr) AliWarning(Form("Failed storing track#%d", index));
273     }
274     
275     // compute total multiplicity
276     if (rsn->GetMultiplicity() <= 0) {
277         AliWarning("Zero Multiplicity in this event");
278         return kFALSE;
279     }
280     
281     return kTRUE;
282 }
283
284 //_____________________________________________________________________________
285 Bool_t AliRsnReader::FillFromAOD(AliRsnEvent *rsn, AliAODEvent *aod, AliMCEvent *mc)
286 {
287 //
288 // Filler from an AOD event.
289 // Stores all tracks (if a filter is defined, it will store
290 // only the ones which survive the cuts).
291 // If a reference MC event is provided, it is used to store
292 // the MC informations for each track (true PDG code, 
293 // GEANT label of mother, PDG code of mother, if any).
294 // When this is used, the 'source' flag of the output
295 // AliRsnEvent object will be set to 'kAOD'.
296 //
297     
298     // retrieve stack (if possible)
299     AliStack *stack = 0x0;
300     if (mc) stack = mc->Stack();
301     
302     // get number of tracks
303     Int_t ntracks = aod->GetNTracks();
304     if (!ntracks) {
305        AliWarning("No tracks in this event");
306        return kFALSE;
307     }
308     
309     // get primary vertex
310     Double_t vertex[3];
311     vertex[0] = aod->GetPrimaryVertex()->GetX();
312     vertex[1] = aod->GetPrimaryVertex()->GetY();
313     vertex[2] = aod->GetPrimaryVertex()->GetZ();
314     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
315     
316     // store tracks from ESD
317     Int_t  index, label, labmum;
318     Bool_t check;
319     AliAODTrack *aodTrack = 0;
320     AliRsnDaughter temp;
321     TObjArrayIter iter(aod->GetTracks());
322     while ( (aodTrack = (AliAODTrack*)iter.Next()) ) {
323         // retrieve index
324         index = aod->GetTracks()->IndexOf(aodTrack);
325         label = aodTrack->GetLabel();
326         if (fRejectFakes && (label < 0)) continue;
327         // copy ESD track data into RsnDaughter
328         // if unsuccessful, this track is skipped
329         check = temp.Adopt(aodTrack);
330         if (!check) continue;
331         // if stack is present, copy MC info
332         if (stack) {
333             TParticle *part = stack->Particle(TMath::Abs(label));
334             if (part) {
335                 temp.InitMCInfo(part);
336                 labmum = part->GetFirstMother();
337                 if (labmum >= 0) {
338                     TParticle *mum = stack->Particle(labmum);
339                     temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
340                 }
341             }
342         }
343         // set index and label and add this object to the output container
344         temp.SetIndex(index);
345         temp.SetLabel(label);
346         AliRsnDaughter *ptr = rsn->AddTrack(temp);
347         // if problems occurred while storin, that pointer is NULL
348         if (!ptr) AliWarning(Form("Failed storing track#%d"));
349     }
350     
351     // compute total multiplicity
352     if (rsn->GetMultiplicity() <= 0) {
353         AliWarning("Zero multiplicity in this event");
354         return kFALSE;
355     }
356     
357     return kTRUE;
358 }
359
360 //_____________________________________________________________________________
361 Bool_t AliRsnReader::FillFromMC(AliRsnEvent *rsn, AliMCEvent *mc)
362 {
363 //
364 // Filler from an ESD event.
365 // Stores all tracks which generate at least one 
366 // TrackReference (a point in a sensitive volume).
367 // In this case, the MC info is stored by default and 
368 // perfect particle identification is the unique available.
369 // When this is used, the 'source' flag of the output
370 // AliRsnEvent object will be set to 'kMC'.
371 //
372     
373     // get number of tracks
374     Int_t ntracks = mc->GetNumberOfTracks();
375     if (!ntracks) {
376        AliWarning("No tracks in this event");
377        return kFALSE;
378     }
379     
380     AliStack *stack = mc->Stack();
381     
382     // get primary vertex
383     TArrayF fvertex(3);
384     Double_t vertex[3];
385     mc->GenEventHeader()->PrimaryVertex(fvertex);
386     vertex[0] = (Double_t)fvertex[0];
387     vertex[1] = (Double_t)fvertex[1];
388     vertex[2] = (Double_t)fvertex[2];
389     rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
390     
391     // store tracks from MC
392     Int_t  index, labmum;
393     Bool_t check;
394     AliRsnDaughter temp;
395     for (index = 0; index < ntracks; index++) {
396         // get and check MC track
397         AliMCParticle *mcTrack = mc->GetTrack(index);
398         // if particle has no track references, it is rejected
399         if (mcTrack->GetNumberOfTrackReferences() <= 0) continue;
400         // try to insert in the RsnDaughter its data
401         check = temp.Adopt(mcTrack);
402         if (!check) continue;
403         labmum = temp.GetMCInfo()->Mother();
404         if (labmum >= 0) {
405             TParticle *mum = stack->Particle(labmum);
406             temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
407         }
408         // if successful, set other data and stores it
409         temp.SetIndex(index);
410         temp.SetLabel(mcTrack->Label());
411         AliRsnDaughter *ptr = rsn->AddTrack(temp);
412         // if problems occurred while storin, that pointer is NULL
413         if (!ptr) AliWarning(Form("Failed storing track#%d", index));
414     }
415     
416     // compute total multiplicity
417     if (rsn->GetMultiplicity() <= 0) {
418         AliWarning("Zero multiplicity in this event");
419         return kFALSE;
420     }
421     
422     return kTRUE;
423 }