Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliESD.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 /* $Id$ */
17
18 //-----------------------------------------------------------------
19 //           Implementation of the ESD class
20 //   This is the class to deal with during the phisical analysis of data
21 //   This class is generated directly by the reconstruction methods
22 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
23 //-----------------------------------------------------------------
24
25 #include "AliESD.h"
26 #include "AliESDfriend.h"
27
28 ClassImp(AliESD)
29
30 //______________________________________________________________________________
31 AliESD::AliESD():
32   fEventNumberInFile(0),
33   fBunchCrossNumber(0),
34   fOrbitNumber(0),
35   fPeriodNumber(0),
36   fRunNumber(0),
37   fTimeStamp(0),
38   fEventType(0),
39   fTriggerMask(0),
40   fTriggerCluster(0),
41   fRecoVersion(0),
42   fMagneticField(0),
43   fZDCN1Energy(0),
44   fZDCP1Energy(0),
45   fZDCN2Energy(0),
46   fZDCP2Energy(0),
47   fZDCEMEnergy(0),
48   fZDCParticipants(0),
49   fT0zVertex(0),
50   fSPDVertex(),
51   fPrimaryVertex(),
52   fSPDMult(),
53   fT0clock(0),
54   fT0timeStart(0),
55   fT0trig(0),
56   fTracks("AliESDtrack",15000),
57   fHLTConfMapTracks("AliESDHLTtrack",25000),
58   fHLTHoughTracks("AliESDHLTtrack",15000),
59   fMuonTracks("AliESDMuonTrack",30),
60   fPmdTracks("AliESDPmdTrack",3000),
61   fTrdTracks("AliESDTrdTrack",300),
62   fV0s("AliESDv0",200),  
63   fCascades("AliESDcascade",20),
64   fKinks("AliESDkink",4000),
65   fCaloClusters("AliESDCaloCluster",10000),
66   fEMCALClusters(0), 
67   fFirstEMCALCluster(-1),
68   fEMCALTriggerPosition(0x0),
69   fEMCALTriggerAmplitudes(0x0),
70   fPHOSClusters(0), 
71   fFirstPHOSCluster(-1),
72   fPHOSTriggerPosition(0x0),
73   fPHOSTriggerAmplitudes(0x0),
74   fESDFMD(0x0),
75   fESDVZERO(0x0),
76   fESDACORDE(0x0),
77   fErrorLogs("AliRawDataErrorLog",5)
78 {
79   // 
80   // Standar constructor
81   //
82
83    for (Int_t i=0; i<3; i++) fT0TOF[i] = 0;
84  for (Int_t i=0; i<24; i++) {
85     fT0time[i] = 0;
86     fT0amplitude[i] = 0;
87   }
88   for (Int_t i=0; i<2; i++) fDiamondXY[i]=0.;
89   for (Int_t i=0; i<3; i++) fDiamondCovXY[i]=0.;
90 }
91
92 //______________________________________________________________________________
93 AliESD::AliESD(const AliESD& esd):
94   TObject(esd),
95   fEventNumberInFile(esd.fEventNumberInFile),
96   fBunchCrossNumber(esd.fBunchCrossNumber),
97   fOrbitNumber(esd.fOrbitNumber),
98   fPeriodNumber(esd.fPeriodNumber),
99   fRunNumber(esd.fRunNumber),
100   fTimeStamp(esd.fTimeStamp),
101   fEventType(esd.fEventType),
102   fTriggerMask(esd.fTriggerMask),
103   fTriggerCluster(esd.fTriggerCluster),
104   fRecoVersion(esd.fRecoVersion),
105   fMagneticField(esd.fMagneticField),
106   fZDCN1Energy(esd.fZDCN1Energy),
107   fZDCP1Energy(esd.fZDCP1Energy),
108   fZDCN2Energy(esd.fZDCN2Energy),
109   fZDCP2Energy(esd.fZDCP2Energy),
110   fZDCEMEnergy(esd.fZDCEMEnergy),
111   fZDCParticipants(esd.fZDCParticipants),
112   fT0zVertex(esd.fT0zVertex),
113   fSPDVertex(esd.fSPDVertex),
114   fPrimaryVertex(esd.fPrimaryVertex),
115   fSPDMult(esd.fSPDMult),
116   fT0clock(esd.fT0clock),
117   fT0timeStart(esd.fT0timeStart),
118   fT0trig(esd.fT0trig),
119   fTracks(*((TClonesArray*)esd.fTracks.Clone())),
120   fHLTConfMapTracks(*((TClonesArray*)esd.fHLTConfMapTracks.Clone())),
121   fHLTHoughTracks(*((TClonesArray*)esd.fHLTHoughTracks.Clone())),
122   fMuonTracks(*((TClonesArray*)esd.fMuonTracks.Clone())),
123   fPmdTracks(*((TClonesArray*)esd.fPmdTracks.Clone())),
124   fTrdTracks(*((TClonesArray*)esd.fTrdTracks.Clone())),
125   fV0s(*((TClonesArray*)esd.fV0s.Clone())),  
126   fCascades(*((TClonesArray*)esd.fCascades.Clone())),
127   fKinks(*((TClonesArray*)esd.fKinks.Clone())),
128   fCaloClusters(*((TClonesArray*)esd.fCaloClusters.Clone())),
129   fEMCALClusters(esd.fEMCALClusters), 
130   fFirstEMCALCluster(esd.fFirstEMCALCluster),
131   fEMCALTriggerPosition(esd. fEMCALTriggerPosition),
132   fEMCALTriggerAmplitudes(esd.fEMCALTriggerAmplitudes),
133   fPHOSClusters(esd.fPHOSClusters), 
134   fFirstPHOSCluster(esd.fFirstPHOSCluster),
135   fPHOSTriggerPosition(esd.fPHOSTriggerPosition),
136   fPHOSTriggerAmplitudes(esd.fPHOSTriggerAmplitudes),
137   fESDFMD(esd.fESDFMD),
138   fESDVZERO(esd.fESDVZERO),
139   fESDACORDE(esd.fESDACORDE),
140   fErrorLogs(*((TClonesArray*)esd.fErrorLogs.Clone()))
141 {
142   // 
143   // copy constructor
144   //
145   for (Int_t i=0; i<3; i++)fT0TOF[i] = esd.fT0TOF[i];
146   for (Int_t i=0; i<24; i++) {
147     fT0time[i] = esd.fT0time[i];
148     fT0amplitude[i] = esd.fT0amplitude[i];
149   }
150   for (Int_t i=0; i<2; i++) fDiamondXY[i]=esd.fDiamondXY[i];
151   for (Int_t i=0; i<3; i++) fDiamondCovXY[i]=esd.fDiamondCovXY[i];
152 }
153
154 //______________________________________________________________________________
155 AliESD::~AliESD()
156 {
157   //
158   // Standard destructor
159   //
160   fTracks.Delete();
161   fHLTConfMapTracks.Delete();
162   fHLTHoughTracks.Delete();
163   fMuonTracks.Delete();
164   fPmdTracks.Delete();
165   fTrdTracks.Delete();
166   fV0s.Delete();
167   fCascades.Delete();
168   fKinks.Delete();
169   fCaloClusters.Delete();
170   delete fESDFMD;
171   delete fESDVZERO;
172   delete fEMCALTriggerPosition;
173   delete fEMCALTriggerAmplitudes;
174   delete fPHOSTriggerPosition;
175   delete fPHOSTriggerAmplitudes;
176   delete fESDACORDE;
177
178   fErrorLogs.Delete();
179
180 }
181
182 //______________________________________________________________________________
183 void AliESD::Reset()
184 {
185   // 
186   // Reset the contents and delete the entries in TClonesArrays
187   //
188
189   fEventNumberInFile=0;
190   fBunchCrossNumber=0;
191   fOrbitNumber=0;
192   fPeriodNumber=0;
193   fRunNumber=0;
194   fTimeStamp = 0;
195   fEventType = 0;
196   fTriggerMask=0;
197   fTriggerCluster=0;
198   fRecoVersion=0;
199   fMagneticField=0;
200   fZDCN1Energy=0;
201   fZDCP1Energy=0;
202   fZDCN2Energy=0;
203   fZDCP2Energy=0;
204   fZDCEMEnergy=0;
205   fZDCParticipants=0;
206   fT0zVertex=0;
207   
208   for (Int_t i=0; i<2; i++) fDiamondXY[i]=0.;
209   for (Int_t i=0; i<3; i++) fDiamondCovXY[i]=0.;
210
211   for (Int_t i=0; i<24; i++) {
212     fT0time[i] = 0;
213     fT0amplitude[i] = 0;
214   }
215   fT0timeStart = 0;
216   fT0clock = 0;
217   for (Int_t i=0; i<3; i++) fT0TOF[i] = 0;
218 //
219   fSPDMult.~AliMultiplicity();
220   new (&fSPDMult) AliMultiplicity();
221   fSPDVertex.~AliESDVertex();
222   new (&fSPDVertex) AliESDVertex();
223   fPrimaryVertex.~AliESDVertex();
224   new (&fPrimaryVertex) AliESDVertex();
225 //
226   fTracks.Delete();
227   fHLTConfMapTracks.Delete();
228   fHLTHoughTracks.Delete();
229   fMuonTracks.Delete();
230   fPmdTracks.Delete();
231   fTrdTracks.Delete();
232   fV0s.Delete();
233   fCascades.Delete();
234   fKinks.Delete();
235   fCaloClusters.Delete();
236 //
237   fEMCALClusters=0; 
238   fFirstEMCALCluster=-1; 
239   fPHOSClusters=0; 
240   fFirstPHOSCluster=-1; 
241 //
242   if (fEMCALTriggerPosition)   fEMCALTriggerPosition  ->Reset();
243   if (fEMCALTriggerAmplitudes) fEMCALTriggerAmplitudes->Reset();
244   if (fPHOSTriggerPosition)    fPHOSTriggerPosition   ->Reset();
245   if (fPHOSTriggerAmplitudes)  fPHOSTriggerAmplitudes ->Reset();
246 //
247   if (fESDFMD) fESDFMD->Clear();
248 //
249   if (fESDVZERO){
250       fESDVZERO->~AliESDVZERO();
251       new (fESDVZERO) AliESDVZERO();
252   } 
253 //
254   if (fESDACORDE){
255       fESDACORDE->~AliESDACORDE();
256       new (fESDACORDE) AliESDACORDE();  
257   }
258 //
259   fErrorLogs.Delete();
260 }
261
262
263 Bool_t  AliESD::RemoveKink(Int_t rm) {
264   // ---------------------------------------------------------
265   // Remove a kink candidate and references to it from ESD,
266   // if this candidate does not come from a reconstructed decay
267   // Not yet implemented...
268   // ---------------------------------------------------------
269   Int_t last=GetNumberOfKinks()-1;
270   if ((rm<0)||(rm>last)) return kFALSE;
271
272   return kTRUE;
273 }
274
275 Bool_t  AliESD::RemoveV0(Int_t rm) {
276   // ---------------------------------------------------------
277   // Remove a V0 candidate and references to it from ESD,
278   // if this candidate does not come from a reconstructed decay
279   // ---------------------------------------------------------
280   Int_t last=GetNumberOfV0s()-1;
281   if ((rm<0)||(rm>last)) return kFALSE;
282
283   AliESDv0 *v0=GetV0(rm);
284   Int_t idxP=v0->GetPindex(), idxN=v0->GetNindex();
285
286   v0=GetV0(last);
287   Int_t lastIdxP=v0->GetPindex(), lastIdxN=v0->GetNindex();
288
289   Int_t used=0;
290
291   // Check if this V0 comes from a reconstructed decay
292   Int_t ncs=GetNumberOfCascades();
293   for (Int_t n=0; n<ncs; n++) {
294     AliESDcascade *cs=GetCascade(n);
295
296     Int_t csIdxP=cs->GetPindex();
297     Int_t csIdxN=cs->GetNindex();
298
299     if (idxP==csIdxP)
300        if (idxN==csIdxN) return kFALSE;
301
302     if (csIdxP==lastIdxP)
303        if (csIdxN==lastIdxN) used++;
304   }
305
306   //Replace the removed V0 with the last V0 
307   TClonesArray &a=fV0s;
308   delete a.RemoveAt(rm);
309
310   if (rm==last) return kTRUE;
311
312   //v0 is pointing to the last V0 candidate... 
313   new (a[rm]) AliESDv0(*v0);
314   delete a.RemoveAt(last);
315
316   if (!used) return kTRUE;
317   
318
319   // Remap the indices of the daughters of reconstructed decays
320   for (Int_t n=0; n<ncs; n++) {
321     AliESDcascade *cs=GetCascade(n);
322
323
324     Int_t csIdxP=cs->GetPindex();
325     Int_t csIdxN=cs->GetNindex();
326
327     if (csIdxP==lastIdxP)
328       if (csIdxN==lastIdxN) {
329          cs->AliESDv0::SetIndex(1,idxP);
330          cs->AliESDv0::SetIndex(0,idxN);
331          used--;
332          if (!used) return kTRUE;
333       }
334   }
335
336   return kTRUE;
337 }
338
339 Bool_t  AliESD::RemoveTrack(Int_t rm) {
340   // ---------------------------------------------------------
341   // Remove a track and references to it from ESD,
342   // if this track does not come from a reconstructed decay
343   // ---------------------------------------------------------
344   Int_t last=GetNumberOfTracks()-1;
345   if ((rm<0)||(rm>last)) return kFALSE;
346
347   Int_t used=0;
348
349   // Check if this track comes from the reconstructed primary vertex
350   if (fPrimaryVertex.GetStatus()) {
351      UShort_t *primIdx=fPrimaryVertex.GetIndices();
352      Int_t n=fPrimaryVertex.GetNIndices();
353      while (n--) {
354        Int_t idx=Int_t(primIdx[n]);
355        if (rm==idx) return kFALSE;
356        if (idx==last) used++; 
357      }
358   }
359   
360   // Check if this track comes from a reconstructed decay
361   Int_t nv0=GetNumberOfV0s();
362   for (Int_t n=0; n<nv0; n++) {
363     AliESDv0 *v0=GetV0(n);
364
365     Int_t idx=v0->GetNindex();
366     if (rm==idx) return kFALSE;
367     if (idx==last) used++;
368
369     idx=v0->GetPindex();
370     if (rm==idx) return kFALSE;
371     if (idx==last) used++;
372   }
373
374   Int_t ncs=GetNumberOfCascades();
375   for (Int_t n=0; n<ncs; n++) {
376     AliESDcascade *cs=GetCascade(n);
377
378     Int_t idx=cs->GetIndex();
379     if (rm==idx) return kFALSE;
380     if (idx==last) used++;
381   }
382
383   Int_t nkn=GetNumberOfKinks();
384   for (Int_t n=0; n<nkn; n++) {
385     AliESDkink *kn=GetKink(n);
386
387     Int_t idx=kn->GetIndex(0);
388     if (rm==idx) return kFALSE;
389     if (idx==last) used++;
390
391     idx=kn->GetIndex(1);
392     if (rm==idx) return kFALSE;
393     if (idx==last) used++;
394   }
395
396
397   //Replace the removed track with the last track 
398   TClonesArray &a=fTracks;
399   delete a.RemoveAt(rm);
400
401   if (rm==last) return kTRUE;
402
403   AliESDtrack *t=GetTrack(last);
404   t->SetID(rm);
405   new (a[rm]) AliESDtrack(*t);
406   delete a.RemoveAt(last);
407
408
409   if (!used) return kTRUE;
410   
411
412   // Remap the indices of the tracks used for the primary vertex reconstruction
413   if (fPrimaryVertex.GetStatus()) {
414      UShort_t *primIdx=fPrimaryVertex.GetIndices();
415      Int_t n=fPrimaryVertex.GetNIndices();
416      while (n--) {
417        Int_t idx=Int_t(primIdx[n]);
418        if (idx==last) {
419           primIdx[n]=Short_t(rm); 
420           used--;
421           if (!used) return kTRUE;
422        }
423      }
424   }
425   
426   // Remap the indices of the daughters of reconstructed decays
427   for (Int_t n=0; n<nv0; n++) {
428     AliESDv0 *v0=GetV0(n);
429     if (v0->GetIndex(0)==last) {
430        v0->SetIndex(0,rm);
431        used--;
432        if (!used) return kTRUE;
433     }
434     if (v0->GetIndex(1)==last) {
435        v0->SetIndex(1,rm);
436        used--;
437        if (!used) return kTRUE;
438     }
439   }
440
441   for (Int_t n=0; n<ncs; n++) {
442     AliESDcascade *cs=GetCascade(n);
443     if (cs->GetIndex()==last) {
444        cs->SetIndex(rm);
445        used--;
446        if (!used) return kTRUE;
447     }
448   }
449
450   for (Int_t n=0; n<nkn; n++) {
451     AliESDkink *kn=GetKink(n);
452     if (kn->GetIndex(0)==last) {
453        kn->SetIndex(rm,0);
454        used--;
455        if (!used) return kTRUE;
456     }
457     if (kn->GetIndex(1)==last) {
458        kn->SetIndex(rm,1);
459        used--;
460        if (!used) return kTRUE;
461     }
462   }
463
464   return kTRUE;
465 }
466
467
468 Bool_t AliESD::Clean(Float_t *cleanPars) {
469   //
470   // Remove the data which are not needed for the physics analysis.
471   //
472   // 1) Cleaning the V0 candidates
473   //    ---------------------------
474   //    If the cosine of the V0 pointing angle "csp" and 
475   //    the DCA between the daughter tracks "dca" does not satisfy 
476   //    the conditions 
477   //
478   //     csp > cleanPars[1] + dca/cleanPars[0]*(1.- cleanPars[1])
479   //
480   //    an attempt to remove this V0 candidate from ESD is made.
481   //
482   //    The V0 candidate gets removed if it does not belong to any 
483   //    recosntructed cascade decay
484   //
485   //    12.11.2007, optimal values: cleanPars[0]=0.5, cleanPars[1]=0.999
486   //
487   // 2) Cleaning the tracks
488   //    ----------------------
489   //    If track's transverse parameter is larger than cleanPars[2]
490   //                       OR
491   //    track's longitudinal parameter is larger than cleanPars[3]
492   //    an attempt to remove this track from ESD is made.
493   //
494   //    The track gets removed if it does not come 
495   //    from a reconstructed decay
496   //
497   Bool_t rc=kFALSE;
498
499   Float_t dcaMax=cleanPars[0];
500   Float_t cspMin=cleanPars[1];
501
502   Int_t nV0s=GetNumberOfV0s();
503   for (Int_t i=nV0s-1; i>=0; i--) {
504     AliESDv0 *v0=GetV0(i);
505
506     Float_t dca=v0->GetDcaV0Daughters();
507     Float_t csp=v0->GetV0CosineOfPointingAngle();
508     Float_t cspcut=cspMin + dca/dcaMax*(1.-cspMin);
509     if (csp > cspcut) continue;
510
511     if (RemoveV0(i)) rc=kTRUE;
512   }
513
514
515   Float_t dmax=cleanPars[2], zmax=cleanPars[3];
516
517   const AliESDVertex *vertex=GetVertex();
518   Bool_t vtxOK=vertex->GetStatus();
519   
520   Int_t nTracks=GetNumberOfTracks();
521   for (Int_t i=nTracks-1; i>=0; i--) {
522     AliESDtrack *track=GetTrack(i);
523     Float_t xy,z; track->GetImpactParameters(xy,z);
524     if ((TMath::Abs(xy) > dmax) || (vtxOK && (TMath::Abs(z) > zmax))) {
525       if (RemoveTrack(i)) rc=kTRUE;
526     }
527   }
528
529   return rc;
530 }
531
532 Int_t AliESD::AddV0(const AliESDv0 *v) {
533   //
534   // Add V0
535   //
536     Int_t idx=fV0s.GetEntriesFast();
537     new(fV0s[idx]) AliESDv0(*v);
538     return idx;
539 }  
540
541 //______________________________________________________________________________
542 void AliESD::Print(Option_t *) const 
543 {
544   //
545   // Print header information of the event
546   //
547   printf("ESD run information\n");
548   printf("Event # in file %d Bunch crossing # %d Orbit # %d Period # %d Run # %d Trigger %lld Magnetic field %f \n",
549          GetEventNumberInFile(),
550          GetBunchCrossNumber(),
551          GetOrbitNumber(),
552          GetPeriodNumber(),
553          GetRunNumber(),
554          GetTriggerMask(),
555          GetMagneticField() );
556     printf("Vertex: (%.4f +- %.4f, %.4f +- %.4f, %.4f +- %.4f) cm\n",
557            fPrimaryVertex.GetXv(), fPrimaryVertex.GetXRes(),
558            fPrimaryVertex.GetYv(), fPrimaryVertex.GetYRes(),
559            fPrimaryVertex.GetZv(), fPrimaryVertex.GetZRes());
560     printf("Mean vertex in RUN: X=%.4f Y=%.4f cm\n",
561            GetDiamondX(),GetDiamondY());
562     printf("SPD Multiplicity. Number of tracklets %d \n",
563            fSPDMult.GetNumberOfTracklets());
564   printf("Event from reconstruction version %d \n",fRecoVersion);
565   printf("Number of tracks: \n");
566   printf("                 charged   %d\n", GetNumberOfTracks());
567   printf("                 hlt CF    %d\n", GetNumberOfHLTConfMapTracks());
568   printf("                 hlt HT    %d\n", GetNumberOfHLTHoughTracks());
569   printf("                 muon      %d\n", GetNumberOfMuonTracks());
570   printf("                 pmd       %d\n", GetNumberOfPmdTracks());
571   printf("                 trd       %d\n", GetNumberOfTrdTracks());
572   printf("                 v0        %d\n", GetNumberOfV0s());
573   printf("                 cascades  %d\n", GetNumberOfCascades());
574   printf("                 kinks     %d\n", GetNumberOfKinks());
575   printf("                 CaloClusters %d\n", GetNumberOfCaloClusters());
576   printf("                 phos      %d\n", GetNumberOfPHOSClusters());
577   printf("                 emcal     %d\n", GetNumberOfEMCALClusters());
578   printf("                 FMD       %s\n", (fESDFMD ? "yes" : "no"));
579   printf("                 VZERO     %s\n", (fESDVZERO ? "yes" : "no"));
580 }
581
582 void AliESD::SetESDfriend(const AliESDfriend *ev) {
583   //
584   // Attaches the complementary info to the ESD
585   //
586   if (!ev) return;
587
588   Int_t ntrk=ev->GetNumberOfTracks();
589
590   for (Int_t i=0; i<ntrk; i++) {
591     const AliESDfriendTrack *f=ev->GetTrack(i);
592     GetTrack(i)->SetFriendTrack(f);
593   }
594 }
595
596 void AliESD::GetESDfriend(AliESDfriend *ev) const {
597   //
598   // Extracts the complementary info from the ESD
599   //
600   if (!ev) return;
601
602   Int_t ntrk=GetNumberOfTracks();
603
604   for (Int_t i=0; i<ntrk; i++) {
605     AliESDtrack *t=GetTrack(i);
606     const AliESDfriendTrack *f=t->GetFriendTrack();
607     ev->AddTrack(f);
608
609     t->ReleaseESDfriendTrack();// Not to have two copies of "friendTrack"
610
611   }
612 }
613
614 void AliESD::SetDiamond(const AliESDVertex *vertex)
615 {
616   //
617   // Set the interaction diamond
618   //  
619     fDiamondXY[0]=vertex->GetXv();
620     fDiamondXY[1]=vertex->GetYv();
621     Double_t cov[6];
622     vertex->GetCovMatrix(cov);
623     fDiamondCovXY[0]=cov[0];
624     fDiamondCovXY[1]=cov[1];
625     fDiamondCovXY[2]=cov[2];
626   }