PWG4/CaloCalib/AliAnalysisTaskEMCALClusterizeFast.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 // --- Root ---
  19 #include <TClonesArray.h>
  20 #include <TGeoManager.h>
  21 #include <TObjArray.h>
  22 #include <TString.h>
  23 #include <TTree.h>
  24 #include "AliAODCaloCluster.h"
  25 #include "AliAODEvent.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliCDBEntry.h"
  28 #include "AliCDBManager.h"
  29 #include "AliCaloCalibPedestal.h"
  30 #include "AliEMCALAfterBurnerUF.h"
  31 #include "AliEMCALCalibData.h"
  32 #include "AliEMCALClusterizerNxN.h"
  33 #include "AliEMCALClusterizerv1.h"
  34 #include "AliEMCALDigit.h"
  35 #include "AliEMCALGeometry.h"
  36 #include "AliEMCALRecParam.h"
  37 #include "AliEMCALRecPoint.h"
  38 #include "AliEMCALRecoUtils.h"
  39 #include "AliESDEvent.h"
  40 #include "AliLog.h"
  41
  42 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  43
  44 ClassImp(AliAnalysisTaskEMCALClusterizeFast)
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
  48   : AliAnalysisTaskSE(name),
  49     fRun(0),
  50     fDigitsArr(0),
  51     fClusterArr(0),
  52     fRecParam(new AliEMCALRecParam),
  53     fClusterizer(0),
  54     fUnfolder(0),
  55     fJustUnfold(kFALSE),
  56     fGeomName("EMCAL_FIRSTYEARV1"),
  57     fGeomMatrixSet(kFALSE),
  58     fLoadGeomMatrices(kFALSE),
  59     fOCDBpath(),
  60     fCalibData(0),
  61     fPedestalData(0),
  62     fOutputAODBranch(0),
  63     fOutputAODBrName(),
  64     fRecoUtils(0),
  65     fLoadCalib(0),
  66     fLoadPed(0)
  67 {
  68   // Constructor
  69
  70   fBranchNames     = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
  71   for(Int_t i = 0; i < 10; ++i)
  72     fGeomMatrix[i] = 0;
  73 }
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
  77 {
  78   // Destructor.
  79
  80   delete fDigitsArr;
  81   delete fClusterizer;
  82   delete fUnfolder;
  83   delete fRecoUtils;
  84 }
  85
  86 //-------------------------------------------------------------------
  87 void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
  88 {
  89   // Create output objects.
  90
  91   if (!fOutputAODBrName.IsNull()) {
  92     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
  93     fOutputAODBranch->SetName(fOutputAODBrName);
  94     AddAODBranch("TClonesArray", &fOutputAODBranch);
  95     AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
  96   }
  97 }
  98
  99 //________________________________________________________________________
 100 void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
 101 {
 102   // Main loop, called for each event
 103
 104   // remove the contents of output list set in the previous event
 105   if (fOutputAODBranch)
 106     fOutputAODBranch->Clear("C");
 107
 108   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 109   AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
 110
 111   if (!esdevent&&!aodevent) {
 112     Error("UserExec","Event not available");
 113     return;
 114   }
 115
 116   LoadBranches();
 117
 118   Init();
 119
 120   if (fJustUnfold) {
 121     AliWarning("Unfolding not implemented");
 122   } else {
 123     if (esdevent)
 124       FillDigitsArray(esdevent);
 125     else
 126       FillDigitsArray(aodevent);
 127     fClusterizer->Digits2Clusters("");
 128     if (esdevent &&fRecoUtils)
 129       fRecoUtils->FindMatches(esdevent,fClusterArr);
 130     if (fOutputAODBranch)
 131       RecPoints2Clusters();
 132   }
 133 }
 134
 135 //________________________________________________________________________
 136 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
 137 {
 138   // Fill digits from cells.
 139
 140   fDigitsArr->Clear("C");
 141   AliAODCaloCells *cells = event->GetEMCALCells();
 142   Int_t ncells = cells->GetNumberOfCells();
 143   if (ncells>fDigitsArr->GetSize())
 144     fDigitsArr->Expand(2*ncells);
 145   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 146     Double_t cellAmplitude=0, cellTime=0;
 147     Short_t cellNumber=0;
 148     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 149       break;
 150     AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
 151     digit->SetId(cellNumber);
 152     digit->SetAmplitude(cellAmplitude);
 153     digit->SetTime(cellTime);
 154     digit->SetTimeR(cellTime);
 155     digit->SetIndexInList(idigit);
 156     digit->SetType(AliEMCALDigit::kHG);
 157     idigit++;
 158   }
 159 }
 160
 161 //________________________________________________________________________
 162 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
 163 {
 164   // Fill digits from cells.
 165
 166   fDigitsArr->Clear("C");
 167   AliESDCaloCells *cells = event->GetEMCALCells();
 168   Int_t ncells = cells->GetNumberOfCells();
 169   if (ncells>fDigitsArr->GetSize())
 170     fDigitsArr->Expand(2*ncells);
 171   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 172     Double_t cellAmplitude=0, cellTime=0;
 173     Short_t cellNumber=0;
 174     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 175       break;
 176     AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
 177     digit->SetId(cellNumber);
 178     digit->SetAmplitude(cellAmplitude);
 179     digit->SetTime(cellTime);
 180     digit->SetTimeR(cellTime);
 181     digit->SetIndexInList(idigit);
 182     digit->SetType(AliEMCALDigit::kHG);
 183     idigit++;
 184   }
 185 }
 186
 187 //________________________________________________________________________________________
 188 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
 189 {
 190   // Cluster energy, global position, cells and their amplitude fractions are restored.
 191
 192   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 193
 194   Int_t Ncls = fClusterArr->GetEntriesFast();
 195   for(Int_t i=0, nout=0; i < Ncls; ++i) {
 196     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 197     Int_t ncells_true = 0;
 198     const Int_t ncells = recpoint->GetMultiplicity();
 199     UShort_t   absIds[ncells];
 200     Double32_t ratios[ncells];
 201     Int_t *dlist = recpoint->GetDigitsList();
 202     Float_t *elist = recpoint->GetEnergiesList();
 203
 204     for (Int_t c = 0; c < ncells; ++c) {
 205       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 206       absIds[ncells_true] = digit->GetId();
 207       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
 208       if (ratios[ncells_true] > 0.001)
 209         ++ncells_true;
 210     }
 211
 212     if (ncells_true < 1) {
 213       AliWarning("Skipping cluster with no cells");
 214       continue;
 215     }
 216
 217     // calculate new cluster position
 218     TVector3 gpos;
 219     Float_t g[3];
 220
 221     recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
 222     recpoint->GetGlobalPosition(gpos);
 223     gpos.GetXYZ(g);
 224
 225     AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(fOutputAODBranch->New(nout++));
 226     clus->SetType(AliVCluster::kEMCALClusterv1);
 227     clus->SetE(recpoint->GetEnergy());
 228     clus->SetPosition(g);
 229     clus->SetNCells(ncells_true);
 230     clus->SetCellsAbsId(absIds);
 231     clus->SetCellsAmplitudeFraction(ratios);
 232     clus->SetDispersion(recpoint->GetDispersion());
 233     clus->SetChi2(-1);                      //not yet implemented
 234     clus->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 235     clus->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 236     Float_t elipAxis[2];
 237     recpoint->GetElipsAxis(elipAxis);
 238     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 239     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 240     clus->SetDistToBadChannel(recpoint->GetDistanceToBadTower());
 241     if (esdevent && fRecoUtils) {
 242       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 243       if(trackIndex >= 0) {
 244         clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
 245         if(DebugLevel() > 1)
 246           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
 247       }
 248     }
 249   }
 250 }
 251
 252 //________________________________________________________________________________________
 253 void AliAnalysisTaskEMCALClusterizeFast::Init()
 254 {
 255   //Select clusterization/unfolding algorithm and set all the needed parameters
 256
 257   AliVEvent * event = InputEvent();
 258   if (!event) {
 259     AliWarning("Event not available!!!");
 260     return;
 261   }
 262
 263   if (event->GetRunNumber()==fRun)
 264     return;
 265   fRun = event->GetRunNumber();
 266
 267   if (fJustUnfold){
 268     // init the unfolding afterburner
 269     delete fUnfolder;
 270     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 271     return;
 272   }
 273
 274   AliCDBManager *cdb = AliCDBManager::Instance();
 275   if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
 276     cdb->SetDefaultStorage(fOCDBpath);
 277   if (fRun!=cdb->GetRun())
 278     cdb->SetRun(fRun);
 279
 280   AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
 281   if (!geometry) {
 282     AliFatal("Geometry not available!!!");
 283     return;
 284   }
 285
 286   if (!fGeomMatrixSet) {
 287     if (fLoadGeomMatrices) {
 288       for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
 289         if(fGeomMatrix[mod]){
 290           if(DebugLevel() > 2)
 291             fGeomMatrix[mod]->Print();
 292           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
 293         }
 294       }
 295     } else {
 296       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
 297         if(event->GetEMCALMatrix(mod)) {
 298           if(DebugLevel() > 2)
 299             event->GetEMCALMatrix(mod)->Print();
 300           geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
 301         }
 302       }
 303     }
 304     fGeomMatrixSet=kTRUE;
 305   }
 306
 307   // setup digit array if needed
 308   if (!fDigitsArr) {
 309     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
 310     fDigitsArr->SetOwner(1);
 311   }
 312
 313   // then setup clusterizer
 314   delete fClusterizer;
 315   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 316     fClusterizer = new AliEMCALClusterizerv1(geometry);
 317   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 318     fClusterizer = new AliEMCALClusterizerNxN(geometry);
 319   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
 320    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
 321    clusterizer->SetNRowDiff(2);
 322    clusterizer->SetNColDiff(2);
 323    fClusterizer = clusterizer;
 324   } else{
 325     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 326   }
 327   fClusterizer->InitParameters(fRecParam);
 328   if (!fCalibData&&fLoadCalib) {
 329     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
 330     if (entry)
 331       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
 332     if (!fCalibData)
 333       AliFatal("Calibration parameters not found in CDB!");
 334   }
 335   if (!fPedestalData&&fLoadPed) {
 336     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
 337     if (entry)
 338       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
 339   }
 340   if (fCalibData) {
 341     fClusterizer->SetInputCalibrated(kFALSE);
 342     fClusterizer->SetCalibrationParameters(fCalibData);
 343     fClusterizer->SetCaloCalibPedestal(fPedestalData);
 344   } else {
 345     fClusterizer->SetInputCalibrated(kTRUE);
 346   }
 347   fClusterizer->SetJustClusters(kTRUE);
 348   fClusterizer->SetDigitsArr(fDigitsArr);
 349   fClusterizer->SetOutput(0);
 350   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 351 }