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