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