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     if (event->GetPHOSCells()) {
 139
 140       //Create the container
 141       if(firstPHOSEvent)
 142       {
 143         fPHOSCells = event->GetPHOSCells() ;// Just recover the pointer of the type ESD or AOD
 144
 145         //Get the arrays and put them in a temporary array
 146         Int_t ncells = fPHOSCells->GetNumberOfCells() ;
 147         Short_t    *cellNumber = new Short_t   [ncells];
 148         Double32_t *amplitude  = new Double32_t[ncells];
 149         Double32_t *time       = new Double32_t[ncells];
 150         for (Int_t icell = 0; icell < ncells; icell++) {
 151           cellNumber[icell] = fPHOSCells->GetCellNumber(icell);
 152           amplitude [icell] = fPHOSCells->GetAmplitude(icell);
 153           time      [icell] = fPHOSCells->GetTime(icell) ;
 154         }
 155
 156         //Now clean the event and enlarge the arrays
 157         fPHOSCells->DeleteContainer();      // Delete the arrays set for the first event we need a larger array
 158         fPHOSCells->SetType(AliVCaloCells::kPHOSCell) ;
 159         fPHOSCells->CreateContainer(fNumberOfPHOSCells) ;
 160
 161         //Put back the contents of the first event to the enlarged array
 162         for (Int_t icell = 0; icell < ncells; icell++) {
 163           fPHOSCells->SetCell(phosPos++, cellNumber[icell], amplitude[icell], time[icell]) ;
 164         }
 165
 166         delete [] cellNumber ;
 167         delete [] amplitude ;
 168         delete [] time ;
 169
 170         firstPHOSEvent=kFALSE;
 171       }
 172       else // Add the rest of the events
 173       {
 174         Int_t ncells = event->GetPHOSCells()->GetNumberOfCells() ;
 175         AliVCaloCells * phosCells = event->GetPHOSCells() ;
 176         for (Int_t icell = 0; icell < ncells; icell++) {
 177           fPHOSCells->SetCell(phosPos++, phosCells->GetCellNumber(icell), phosCells->GetAmplitude(icell), phosCells->GetTime(icell)) ;
 178         }
 179       }//not first event
 180     }
 181     if (event->GetEMCALCells()) {
 182
 183       //Create the container
 184       if(firstEMCALEvent)
 185       {
 186         fEMCALCells = event->GetEMCALCells() ; // Just recover the pointer of the type ESD or AOD
 187
 188         //Get the arrays and put them in a temporary array
 189         Int_t ncells = fEMCALCells->GetNumberOfCells() ;
 190         Short_t    *cellNumber = new Short_t   [ncells];
 191         Double32_t *amplitude  = new Double32_t[ncells];
 192         Double32_t *time       = new Double32_t[ncells];
 193         for (Int_t icell = 0; icell < ncells; icell++) {
 194           cellNumber[icell] = fEMCALCells->GetCellNumber(icell);
 195           amplitude [icell] = fEMCALCells->GetAmplitude(icell);
 196           time      [icell] = fEMCALCells->GetTime(icell) ;
 197         }
 198
 199         //Now clean the event and enlarge the arrays
 200         fEMCALCells->DeleteContainer();        // Delete the arrays set for the first event we need a larger array
 201         fEMCALCells->SetType(AliVCaloCells::kEMCALCell) ;
 202         fEMCALCells->CreateContainer(fNumberOfEMCALCells) ;
 203
 204         //Put back the contents of the first event to the enlarged array
 205         for (Int_t icell = 0; icell < ncells; icell++) {
 206           fEMCALCells->SetCell(emcalPos++, cellNumber[icell], amplitude[icell], time[icell]) ;
 207         }
 208
 209         delete [] cellNumber ;
 210         delete [] amplitude ;
 211         delete [] time ;
 212
 213         firstEMCALEvent=kFALSE;
 214       }
 215       else // Add the rest of the events
 216       {
 217         Int_t ncells = event->GetEMCALCells()->GetNumberOfCells() ;
 218         AliVCaloCells * emcalCells = event->GetEMCALCells() ;
 219         for (Int_t icell = 0; icell < ncells; icell++) {
 220           fEMCALCells->SetCell(emcalPos++, emcalCells->GetCellNumber(icell), emcalCells->GetAmplitude(icell), emcalCells->GetTime(icell)) ;
 221         }
 222       }//not first event
 223     }
 224   }
 225 }
 226
 227 AliVParticle* AliMixedEvent::GetTrack(Int_t i) const
 228 {
 229     // Return track # i
 230     Int_t iEv  = TMath::BinarySearch(fNEvents, fNTracksCumul, i);
 231     while((iEv < (fNEvents - 1)) && (fNTracksCumul[iEv] == fNTracksCumul[iEv+1])) {iEv++;}
 232
 233     Int_t irel = i - fNTracksCumul[iEv];
 234     AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
 235     return (evt->GetTrack(irel));
 236 }
 237
 238 AliVCluster* AliMixedEvent::GetCaloCluster(Int_t i) const
 239 {
 240     // Return calo cluster # i
 241   Int_t iEv  = TMath::BinarySearch(fNEvents, fNCaloClustersCumul, i);
 242   while((iEv < (fNEvents - 1)) && (fNCaloClustersCumul[iEv] == fNCaloClustersCumul[iEv+1])) {iEv++;}
 243
 244   Int_t irel = i - fNCaloClustersCumul[iEv];
 245   AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
 246   return (evt->GetCaloCluster(irel));
 247 }
 248
 249 const AliVVertex* AliMixedEvent::GetEventVertex(Int_t i) const
 250 {
 251     // Return vertex of track # i
 252     Int_t iEv  = TMath::BinarySearch(fNEvents, fNTracksCumul, i);
 253     while((iEv < (fNEvents - 1)) && (fNTracksCumul[iEv] == fNTracksCumul[iEv+1])) {iEv++;}
 254     AliVEvent* evt = (AliVEvent*) (fEventList.At(iEv));
 255     return (evt->GetPrimaryVertex());
 256 }
 257
 258 const AliVVertex* AliMixedEvent::GetVertexOfEvent(Int_t i) const
 259 {
 260     // Return vertex of event # i
 261   if (i > fNEvents)
 262     AliFatal(Form("%d events in buffer, event %d requested", fNEvents, i)) ;
 263   AliVEvent* evt = (AliVEvent*) (fEventList.At(i));
 264   return (evt->GetPrimaryVertex());
 265 }
 266
 267 void AliMixedEvent::Reset()
 268 {
 269     // Reset the event
 270   fEventList.Clear();
 271   fNEvents = 0;
 272   fNumberOfTracks = 0;
 273   fNumberOfCaloClusters = 0;
 274   fNumberOfPHOSCells = 0;
 275   fNumberOfEMCALCells = 0;
 276   if (fNTracksCumul) {
 277     delete[]  fNTracksCumul;
 278     fNTracksCumul = 0;
 279   }
 280   if (fNCaloClustersCumul) {
 281     delete[]  fNCaloClustersCumul;
 282     fNCaloClustersCumul = 0;
 283   }
 284   if (fNPHOSCellsCumul) {
 285     delete[]  fNPHOSCellsCumul;
 286     fNPHOSCellsCumul = 0;
 287   }
 288   if (fNEMCALCellsCumul) {
 289     delete[]  fNEMCALCellsCumul;
 290     fNEMCALCellsCumul = 0;
 291   }
 292   if (fPHOSCells) {
 293      fPHOSCells->DeleteContainer();
 294   }
 295   if (fEMCALCells) {
 296     fEMCALCells->DeleteContainer();
 297   }
 298 }
 299
 300 Int_t AliMixedEvent::EventIndex(Int_t itrack) const
 301 {
 302   // Return the event index for track #itrack
 303   return  TMath::BinarySearch(fNEvents, fNTracksCumul, itrack);
 304 }
 305
 306 Int_t AliMixedEvent::EventIndexForCaloCluster(Int_t icluster) const
 307 {
 308     // Return the event index for track #itrack
 309   return  TMath::BinarySearch(fNEvents, fNCaloClustersCumul, icluster);
 310 }
 311
 312 Int_t AliMixedEvent::EventIndexForPHOSCell(Int_t icell) const
 313 {
 314     // Return the event index for track #itrack
 315   return  TMath::BinarySearch(fNEvents, fNPHOSCellsCumul, icell);
 316 }
 317
 318 Int_t AliMixedEvent::EventIndexForEMCALCell(Int_t icell) const
 319 {
 320     // Return the event index for track #itrack
 321   return  TMath::BinarySearch(fNEvents, fNEMCALCellsCumul, icell);
 322 }
 323
 324 Bool_t AliMixedEvent::ComputeVtx(TObjArray *vertices, Double_t *pos,Double_t *sig,Int_t *nContributors){
 325 //
 326 // Calculate the mean vertex psoitions from events in the buffer
 327
 328     Int_t nentries = vertices->GetEntriesFast();
 329     Double_t sum[3]={0.,0.,0.};
 330     Double_t sumsigma[6]={0.,0.,0.,0.,0.,0.};
 331
 332
 333     for(Int_t ivtx = 0; ivtx < nentries; ivtx++){
 334         AliVVertex *vtx=(AliVVertex*)vertices->UncheckedAt(ivtx);
 335         Double_t covariance[6];
 336         vtx->GetCovarianceMatrix(covariance);
 337         Double_t vtxPos[3];
 338         vtx->GetXYZ(vtxPos);
 339         if(TMath::Abs(covariance[0])<1.e-13) {
 340         return kFALSE;
 341         }else{
 342         sum[0]+=vtxPos[0]*(1./covariance[0]);
 343         sumsigma[0]+=(1./covariance[0]);
 344         }
 345         if(TMath::Abs(covariance[2])<1.e-13) {
 346         return kFALSE;
 347         }else{
 348         sum[1]+=vtxPos[1]*(1./covariance[2]);
 349         sumsigma[2]+=(1./covariance[2]);
 350         }
 351         if(TMath::Abs(covariance[5])<1.e-13) {
 352         return kFALSE;
 353         }else{
 354         sum[2]+=vtxPos[2]*(1./covariance[5]);
 355         sumsigma[5]+=(1./covariance[5]);
 356         }
 357         if(TMath::Abs(covariance[1])<1.e-13) {
 358          sumsigma[1]+=0.;
 359         }else{
 360         sumsigma[1]+=(1./covariance[1]);
 361         }
 362         if(TMath::Abs(covariance[3])<1.e-13) {
 363         sumsigma[3]+=0.;
 364         }else{
 365         sumsigma[3]+=(1./covariance[3]);
 366         }
 367         if(TMath::Abs(covariance[4])<1.e-13) {
 368         sumsigma[4]+=0.;
 369         }else{
 370         sumsigma[4]+=(1./covariance[4]);
 371         }
 372
 373      nContributors[0]=nContributors[0]+vtx->GetNContributors();
 374     }
 375
 376     for(Int_t i=0;i<3;i++){
 377         if(TMath::Abs(sumsigma[i])<1.e-13) continue;
 378         pos[i]=sum[i]/sumsigma[i];
 379     }
 380     for(Int_t i2=0;i2<6;i2++){
 381         if(TMath::Abs(sumsigma[i2])<1.e-13) {sig[i2]=0.; continue;}
 382         sig[i2]=1./sumsigma[i2];
 383     }
 384     return kTRUE;
 385 }
 386
 387
 388 Double_t AliMixedEvent::GetMagneticField() const
 389 {
 390     // Return magnetic field of the first event in the list
 391     if (fEventList.GetEntries() == 0) return -999.;
 392
 393     AliVEvent* evt = (AliVEvent*) (fEventList.At(0));
 394     return evt->GetMagneticField();
 395 }