Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliMixedEvent.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 /* $Id$ */
17
18
19 //-------------------------------------------------------------------------
20 //                          Class AliMixedEvent
21 // VEvent which is the container of several VEvents 
22 // Use Case: Event Mixing     
23 // Origin: Andreas Morsch, CERN, Andreas.Morsch@cern.ch 
24 //-------------------------------------------------------------------------
25
26
27 #include "AliMixedEvent.h"
28 #include "AliExternalTrackParam.h"
29 #include "TVector3.h"
30 #include "AliVVertex.h"
31 #include <TMath.h>
32 #include <TMatrix.h>
33 #include <TMatrixD.h>
34 #include "AliLog.h"
35 #include "AliVCaloCells.h"
36
37
38 ClassImp(AliMixedEvent)
39
40
41 AliMixedEvent::AliMixedEvent() :
42   AliVEvent(),
43   fEventList(),
44   fNEvents(0),       
45   fNumberOfTracks(0),
46   fNumberOfCaloClusters(0), 
47   fNumberOfPHOSCells(0), 
48   fNumberOfEMCALCells(0),
49   fNTracksCumul(0),
50   fNCaloClustersCumul(0),
51   fNPHOSCellsCumul(0), 
52   fNEMCALCellsCumul(0), 
53   fPHOSCells(NULL), 
54   fEMCALCells(NULL), 
55   fMeanVertex(0)
56 {
57     // Default constructor
58 }
59
60 AliMixedEvent::AliMixedEvent(const AliMixedEvent& Evnt) :
61   AliVEvent(Evnt),
62   fEventList(),
63   fNEvents(0),
64   fNumberOfTracks(0),
65   fNumberOfCaloClusters(0), 
66   fNumberOfPHOSCells(0), 
67   fNumberOfEMCALCells(0),
68   fNTracksCumul(0),
69   fNCaloClustersCumul(0),
70   fNPHOSCellsCumul(0), 
71   fNEMCALCellsCumul(0), 
72   fPHOSCells(NULL), 
73   fEMCALCells(NULL), 
74   fMeanVertex(0)
75 { } // Copy constructor
76
77 AliMixedEvent& AliMixedEvent::operator=(const AliMixedEvent& vEvnt)
78 { if (this!=&vEvnt) { 
79     AliVEvent::operator=(vEvnt); 
80   }
81   
82   return *this; 
83 }
84
85 AliMixedEvent::~AliMixedEvent() 
86 {
87     // dtor
88   Reset();
89   delete fPHOSCells ; 
90   delete fEMCALCells ; 
91
92
93
94 void AliMixedEvent::AddEvent(AliVEvent* evt)
95 {
96     // Add a new event to the list
97     fEventList.AddLast(evt);
98 }
99
100
101 void AliMixedEvent::Init()
102 {
103     // Initialize meta information
104   fNEvents = fEventList.GetEntries();
105   fNTracksCumul = new Int_t[fNEvents];
106   fNumberOfTracks = 0;
107   fNCaloClustersCumul = new Int_t[fNEvents];
108   fNumberOfCaloClusters = 0;
109   fNumberOfPHOSCells    = 0;  
110   fNumberOfEMCALCells   = 0; 
111   fNPHOSCellsCumul  = new Int_t[fNEvents];
112   fNEMCALCellsCumul = new Int_t[fNEvents];
113
114   TIter next(&fEventList);
115   AliVEvent* event;
116   Int_t iev = 0;
117     
118   while((event = (AliVEvent*)next())) {
119     fNTracksCumul[iev] = fNumberOfTracks;
120     fNumberOfTracks += (event->GetNumberOfTracks());
121     fNCaloClustersCumul[iev] = fNumberOfCaloClusters;
122     fNumberOfCaloClusters += event->GetNumberOfCaloClusters(); 
123     fNPHOSCellsCumul[iev] = fNumberOfPHOSCells;
124     if (event->GetPHOSCells()) 
125       fNumberOfPHOSCells += event->GetPHOSCells()->GetNumberOfCells(); 
126     fNEMCALCellsCumul[iev] = fNumberOfEMCALCells;
127     if (event->GetEMCALCells()) 
128       fNumberOfEMCALCells += event->GetEMCALCells()->GetNumberOfCells(); 
129     iev++ ;  
130   }
131
132   next.Reset() ; 
133   Short_t phosPos = 0, emcalPos = 0; 
134   Int_t firstPHOSEvent  = kTRUE;
135   Int_t firstEMCALEvent = kTRUE;
136   
137   while((event = (AliVEvent*)next())) {
138     AliVCaloCells * phosCells = event->GetPHOSCells() ; 
139     if (phosCells) {
140       
141       //Create the container
142       if(firstPHOSEvent)
143       {
144         if(!fPHOSCells) fPHOSCells = phosCells->CopyCaloCells(kFALSE) ;// Just recover the first event type:  ESD/AOD
145         else fPHOSCells->DeleteContainer(); //delete the previous container 
146         //Now create a new container with the adequate size
147         fPHOSCells->SetType(AliVCaloCells::kPHOSCell) ; 
148         fPHOSCells->CreateContainer(fNumberOfPHOSCells) ;
149         firstPHOSEvent=kFALSE;
150
151       }//First event
152
153       Int_t ncells = event->GetPHOSCells()->GetNumberOfCells() ;
154       for (Int_t icell = 0; icell < ncells; icell++) {
155           fPHOSCells->SetCell(phosPos++, phosCells->GetCellNumber(icell), phosCells->GetAmplitude(icell), phosCells->GetTime(icell)) ; 
156       }
157      
158     }// phos cells
159     
160     AliVCaloCells * emcalCells = event->GetEMCALCells() ; 
161     if (emcalCells) {
162       
163       //Create the container
164       if(firstEMCALEvent)
165       {
166         if(!fEMCALCells)fEMCALCells = emcalCells->CopyCaloCells(kFALSE) ; // Just recover the first event type:  ESD/AOD
167         else fEMCALCells->DeleteContainer();       // delete the previous container
168         //Now create a new container with the adequate size
169         fEMCALCells->SetType(AliVCaloCells::kEMCALCell) ; 
170         fEMCALCells->CreateContainer(fNumberOfEMCALCells) ;
171         firstEMCALEvent=kFALSE;
172       }//First event
173       
174       Int_t ncells = emcalCells->GetNumberOfCells() ;
175       for (Int_t icell = 0; icell < ncells; icell++) {
176           fEMCALCells->SetCell(emcalPos++, emcalCells->GetCellNumber(icell), emcalCells->GetAmplitude(icell), emcalCells->GetTime(icell)) ; 
177       }
178     }//EMCAL cells
179   }//while event
180   
181 }
182
183 AliVParticle* AliMixedEvent::GetTrack(Int_t i) const
184 {
185     // Return track # i
186     Int_t iEv  = TMath::BinarySearch(fNEvents, fNTracksCumul, i);
187     while((iEv < (fNEvents - 1)) && (fNTracksCumul[iEv] == fNTracksCumul[iEv+1])) {iEv++;}
188
189     Int_t irel = i - fNTracksCumul[iEv];
190     AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
191     return (evt->GetTrack(irel));
192 }
193
194 AliVCluster* AliMixedEvent::GetCaloCluster(Int_t i) const
195 {
196     // Return calo cluster # i
197   Int_t iEv  = TMath::BinarySearch(fNEvents, fNCaloClustersCumul, i);
198   while((iEv < (fNEvents - 1)) && (fNCaloClustersCumul[iEv] == fNCaloClustersCumul[iEv+1])) {iEv++;}
199   
200   Int_t irel = i - fNCaloClustersCumul[iEv];
201   AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
202   return (evt->GetCaloCluster(irel));
203 }
204
205 const AliVVertex* AliMixedEvent::GetEventVertex(Int_t i) const
206 {
207     // Return vertex of track # i
208     Int_t iEv  = TMath::BinarySearch(fNEvents, fNTracksCumul, i);
209     while((iEv < (fNEvents - 1)) && (fNTracksCumul[iEv] == fNTracksCumul[iEv+1])) {iEv++;}
210     AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
211     return (evt->GetPrimaryVertex());
212 }
213
214 const AliVVertex* AliMixedEvent::GetVertexOfEvent(Int_t i) const
215 {
216     // Return vertex of event # i
217   if (i > fNEvents)
218     AliFatal(Form("%d events in buffer, event %d requested", fNEvents, i)) ;  
219   AliVEvent* evt = (AliVEvent*) (fEventList.At(i));
220   return (evt->GetPrimaryVertex());
221 }
222
223 void AliMixedEvent::Reset()
224 {
225     // Reset the event
226   fEventList.Clear();
227   fNEvents = 0;
228   fNumberOfTracks = 0;
229   fNumberOfCaloClusters = 0;
230   fNumberOfPHOSCells = 0;
231   fNumberOfEMCALCells = 0;
232   if (fNTracksCumul) {
233     delete[]  fNTracksCumul;
234     fNTracksCumul = 0;
235   }
236   if (fNCaloClustersCumul) {
237     delete[]  fNCaloClustersCumul;
238     fNCaloClustersCumul = 0;
239   }
240   if (fNPHOSCellsCumul) {
241     delete[]  fNPHOSCellsCumul;
242     fNPHOSCellsCumul = 0;
243   }
244   if (fNEMCALCellsCumul) {
245     delete[]  fNEMCALCellsCumul;
246     fNEMCALCellsCumul = 0;
247   }
248   
249   if (fPHOSCells) {      
250     fPHOSCells->DeleteContainer();       
251   }      
252   if (fEMCALCells) {     
253     fEMCALCells->DeleteContainer();      
254   }
255   
256 }
257
258 Int_t AliMixedEvent::EventIndex(Int_t itrack) const
259 {
260   // Return the event index for track #itrack
261   return  TMath::BinarySearch(fNEvents, fNTracksCumul, itrack);
262 }
263
264 Int_t AliMixedEvent::EventIndexForCaloCluster(Int_t icluster) const
265 {
266     // Return the event index for track #itrack
267   return  TMath::BinarySearch(fNEvents, fNCaloClustersCumul, icluster);
268 }
269
270 Int_t AliMixedEvent::EventIndexForPHOSCell(Int_t icell) const
271 {
272     // Return the event index for track #itrack
273   return  TMath::BinarySearch(fNEvents, fNPHOSCellsCumul, icell);
274 }
275
276 Int_t AliMixedEvent::EventIndexForEMCALCell(Int_t icell) const
277 {
278     // Return the event index for track #itrack
279   return  TMath::BinarySearch(fNEvents, fNEMCALCellsCumul, icell);
280 }
281
282 Bool_t AliMixedEvent::ComputeVtx(TObjArray *vertices, Double_t *pos,Double_t *sig,Int_t *nContributors){
283 //
284 // Calculate the mean vertex psoitions from events in the buffer
285  
286     Int_t nentries = vertices->GetEntriesFast();
287     Double_t sum[3]={0.,0.,0.};
288     Double_t sumsigma[6]={0.,0.,0.,0.,0.,0.};
289
290     
291     for(Int_t ivtx = 0; ivtx < nentries; ivtx++){
292         AliVVertex *vtx=(AliVVertex*)vertices->UncheckedAt(ivtx);
293         Double_t covariance[6];
294         vtx->GetCovarianceMatrix(covariance);
295         Double_t vtxPos[3];
296         vtx->GetXYZ(vtxPos);
297         if(TMath::Abs(covariance[0])<1.e-13) {
298         return kFALSE;
299         }else{
300         sum[0]+=vtxPos[0]*(1./covariance[0]);
301         sumsigma[0]+=(1./covariance[0]);
302         }
303         if(TMath::Abs(covariance[2])<1.e-13) {
304         return kFALSE;
305         }else{
306         sum[1]+=vtxPos[1]*(1./covariance[2]);
307         sumsigma[2]+=(1./covariance[2]);
308         }
309         if(TMath::Abs(covariance[5])<1.e-13) {
310         return kFALSE;
311         }else{
312         sum[2]+=vtxPos[2]*(1./covariance[5]);
313         sumsigma[5]+=(1./covariance[5]);
314         }
315         if(TMath::Abs(covariance[1])<1.e-13) {
316          sumsigma[1]+=0.;
317         }else{
318         sumsigma[1]+=(1./covariance[1]);
319         }
320         if(TMath::Abs(covariance[3])<1.e-13) {
321         sumsigma[3]+=0.;
322         }else{
323         sumsigma[3]+=(1./covariance[3]);
324         }
325         if(TMath::Abs(covariance[4])<1.e-13) {
326         sumsigma[4]+=0.;
327         }else{
328         sumsigma[4]+=(1./covariance[4]);
329         }
330
331      nContributors[0]=nContributors[0]+vtx->GetNContributors();
332     }
333     
334     for(Int_t i=0;i<3;i++){
335         if(TMath::Abs(sumsigma[i])<1.e-13) continue;
336         pos[i]=sum[i]/sumsigma[i];
337     }
338     for(Int_t i2=0;i2<6;i2++){
339         if(TMath::Abs(sumsigma[i2])<1.e-13) {sig[i2]=0.; continue;}
340         sig[i2]=1./sumsigma[i2];
341     }
342     return kTRUE;
343 }
344
345
346 Double_t AliMixedEvent::GetMagneticField() const
347 {
348     // Return magnetic field of the first event in the list
349     if (fEventList.GetEntries() == 0) return -999.;
350     
351     AliVEvent* evt = (AliVEvent*) (fEventList.At(0));
352     return evt->GetMagneticField();
353 }