EMCAL/AliEMCALAfterBurnerUF.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 // After-burner for the EMCAL cluster unfolding algorithm
  16 //
  17 // Input: TObjArray  *clusArray -- array of AliVClusters;
  18 //        AliVCaloCells  *cellsEMCAL -- EMCAL cells.
  19 //
  20 // Output is appended to clusArray, the original (unfolded or not) clusters
  21 // are deleted or moved to another position in clusArray.
  22 //
  23 // If you want to use particular geometry, you must initialize it _before_
  24 // creating AliEMCALAfterBurnerUF instance. Add this or similar line to the
  25 // initialization section:
  26 //
  27 //    AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
  28 //
  29 // gGeoManager must be initialized for this code to work! Do it yourself or
  30 // provide geometry.root file in the current directory so that
  31 // AliEMCALAfterBurnerUF will take it by itself.
  32 // How to use:
  33 //
  34 //   // add this lines to the initialization section of your analysis
  35 //   AliEMCALAfterBurnerUF *abuf = new AliEMCALAfterBurnerUF();
  36 //   TObjArray *clusArray = new TObjArray(100);
  37 //
  38 //
  39 //   AliVEvent *event = InputEvent();
  40 //   AliVCaloCells *cellsEMCAL = event->GetEMCALCells();
  41 //
  42 //   for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
  43 //   {
  44 //     AliVCluster *clus = event->GetCaloCluster(i);
  45 //
  46 //     clusArray->Add(clus->Clone());   // NOTE _CLONE_ in this line
  47 //   }
  48 //
  49 //   abuf->UnfoldClusters(clusArray, cellsEMCAL);
  50 //
  51 //   // do an analysis with clusArray
  52 //   // ....
  53 //
  54 //   // prevent memory leak
  55 //   clusArray->Delete();
  56 //
  57 //----
  58 //  Author: Olga Driga (SUBATECH)
  59
  60 // --- ROOT system ---
  61 #include <TObjArray.h>
  62 #include <TClonesArray.h>
  63 #include <TGeoManager.h>
  64
  65 // --- Standard library --
  66
  67 // --- AliRoot header files ---
  68 #include <AliEMCALAfterBurnerUF.h>
  69 #include <AliEMCALGeometry.h>
  70 #include <AliEMCALUnfolding.h>
  71 #include <AliAODCaloCluster.h>
  72 #include <AliVCaloCells.h>
  73 #include <AliEMCALRecPoint.h>
  74 #include <AliEMCALDigit.h>
  75
  76 ClassImp(AliEMCALAfterBurnerUF)
  77
  78 //------------------------------------------------------------------------
  79 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
  80   fGeom(NULL),
  81   fLogWeight(4.5),      // correct?
  82   fECALocMaxCut(0.03),  // value suggested by Adam
  83   fMinECut(0.01),
  84   fRecPoints(NULL),
  85   fDigitsArr(NULL),
  86   fClusterUnfolding(NULL)
  87 {
  88   // Use this constructor, if unsure
  89
  90   Init();
  91 }
  92
  93 //------------------------------------------------------------------------
  94 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
  95   fGeom(NULL),
  96   fLogWeight(logWeight),
  97   fECALocMaxCut(ecaLocMaxCut),
  98   fMinECut(minECut),
  99   fRecPoints(NULL),
 100   fDigitsArr(NULL),
 101   fClusterUnfolding(NULL)
 102 {
 103   // This constructor allows to set parameters
 104   // Recommended values:
 105   //   Float_t logWeight = 4.5, ECALocMaxCut = 0.03
 106
 107   Init();
 108 }
 109
 110 //------------------------------------------------------------------------
 111 void AliEMCALAfterBurnerUF::Init()
 112 {
 113   // After-burner initialization
 114   // Imports geometry.root (if required), creates unfolding class instance
 115   //
 116   // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
 117   //       for this to work, the OCDB geometry must be imported instead
 118
 119   if (!gGeoManager)
 120     Warning("AliEMCALAfterBurnerUF::Init","GeoManager not found, please import the geometry.root file or pass to the geometry the misalignment matrices");
 121 //    TGeoManager::Import("geometry.root");
 122
 123   // required for global cluster position recalculation
 124   if (!gGeoManager)
 125     Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
 126
 127   // initialize geometry, if not yet initialized
 128   if (!AliEMCALGeometry::GetInstance()) {
 129     Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
 130     AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 131   }
 132
 133   // AliEMCALRecPoint is using exactly this call
 134   fGeom = AliEMCALGeometry::GetInstance();
 135   if (!fGeom)
 136     Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
 137
 138   fClusterUnfolding = new AliEMCALUnfolding(fGeom);
 139   fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
 140   fClusterUnfolding->SetThreshold(fMinECut);
 141
 142   // clusters --> recPoints, cells --> digits and back ;)
 143   fRecPoints = new TObjArray(100);
 144   fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
 145 }
 146
 147 //------------------------------------------------------------------------
 148 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
 149 {
 150   //
 151   // destructor
 152   //
 153
 154   if (fClusterUnfolding) delete fClusterUnfolding;
 155
 156   if (fRecPoints) {
 157     fRecPoints->Delete();
 158     delete fRecPoints;
 159   }
 160   if (fDigitsArr) {
 161     fDigitsArr->Clear("C");
 162     delete fDigitsArr;
 163   }
 164 }
 165
 166 //------------------------------------------------------------------------
 167 void AliEMCALAfterBurnerUF::Clear()
 168 {
 169   //Clean the arrays
 170
 171   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
 172   if (fDigitsArr) fDigitsArr->Clear("C");
 173
 174 }
 175 //------------------------------------------------------------------------
 176 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
 177 {
 178   // Restore clusters from recPoints
 179   // Cluster energy, global position, cells and their amplitude fractions are restored
 180   //
 181   // TODO: restore time and other parameters
 182
 183   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
 184   {
 185     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
 186
 187     const Int_t ncells = recPoint->GetMultiplicity();
 188     Int_t ncellsTrue = 0;
 189
 190     // cells and their amplitude fractions
 191     UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
 192     Double32_t ratios[ncells];
 193
 194     for (Int_t c = 0; c < ncells; c++) {
 195       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
 196
 197       absIds[ncellsTrue] = digit->GetId();
 198       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 199
 200       if (ratios[ncellsTrue] > 0.001) ncellsTrue++;
 201     }
 202
 203     if (ncellsTrue < 1) {
 204       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
 205       continue;
 206     }
 207
 208     TVector3 gpos;
 209     Float_t g[3];
 210
 211     // calculate new cluster position
 212     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
 213     recPoint->GetGlobalPosition(gpos);
 214     gpos.GetXYZ(g);
 215
 216     // create a new cluster
 217     AliAODCaloCluster *clus = new AliAODCaloCluster();
 218     clus->SetType(AliVCluster::kEMCALClusterv1);
 219     clus->SetE(recPoint->GetEnergy());
 220     clus->SetPosition(g);
 221     clus->SetNCells(ncellsTrue);
 222     clus->SetCellsAbsId(absIds);
 223     clus->SetCellsAmplitudeFraction(ratios);
 224     // TODO: time not stored
 225     // TODO: some other properties not stored
 226
 227     clusArray->Add(clus);
 228   } // recPoints loop
 229
 230 }
 231
 232 //------------------------------------------------------------------------
 233 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
 234 {
 235   // Unfolds clusters.
 236   //
 237   // Input: TObjArray of clusters, EMCAL cells.
 238   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
 239
 240   Int_t ndigits = 0;
 241
 242   Int_t nclus = clusArray->GetEntriesFast();
 243
 244   /* Fill recPoints with digits
 245   */
 246   for (Int_t i = 0; i < nclus; i++)
 247   {
 248     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
 249     if (!clus->IsEMCAL()) continue;
 250
 251     // new recPoint
 252     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
 253     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 254     fRecPoints->Add(recPoint);
 255
 256     // fill digits
 257     for (Int_t c = 0; c < clus->GetNCells(); c++) {
 258       Int_t absId = clus->GetCellAbsId(c);
 259       Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
 260       Double_t time = cellsEMCAL->GetCellTime(absId);
 261
 262       // NOTE: it is easy to include cells recalibration here:
 263       // amp *= factor;
 264
 265       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
 266
 267       digit->SetId(absId);
 268       digit->SetAmplitude(amp);
 269       digit->SetTime(time);
 270       digit->SetTimeR(time);
 271       digit->SetIndexInList(ndigits);
 272
 273       recPoint->AddDigit(*digit, amp, kFALSE);
 274
 275       ndigits++;
 276     }
 277
 278     // this cluster will be substituted with the result of unfolding
 279     clusArray->RemoveAt(i);
 280     delete clus;
 281   } // cluster loop
 282
 283
 284   /* Peform unfolding
 285   */
 286   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
 287   fClusterUnfolding->MakeUnfolding();
 288
 289   /* Restore clusters from recPoints.
 290   */
 291   RecPoints2Clusters(clusArray);
 292
 293   // clean up
 294   fRecPoints->Delete(); // do not Clear(), it leaks, why?
 295   fDigitsArr->Clear("C");
 296
 297   clusArray->Compress();
 298
 299 }