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   if (fClusterUnfolding) delete fClusterUnfolding;
 151
 152   if (fRecPoints) {
 153     fRecPoints->Delete();
 154     delete fRecPoints;
 155   }
 156   if (fDigitsArr) {
 157     fDigitsArr->Clear("C");
 158     delete fDigitsArr;
 159   }
 160 }
 161
 162 //------------------------------------------------------------------------
 163 void AliEMCALAfterBurnerUF::Clear()
 164 {
 165   //Clean the arrays
 166
 167   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
 168   if (fDigitsArr) fDigitsArr->Clear("C");
 169
 170 }
 171 //------------------------------------------------------------------------
 172 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
 173 {
 174   // Restore clusters from recPoints
 175   // Cluster energy, global position, cells and their amplitude fractions are restored
 176   //
 177   // TODO: restore time and other parameters
 178
 179   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
 180   {
 181     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
 182
 183     const Int_t ncells = recPoint->GetMultiplicity();
 184     Int_t ncells_true = 0;
 185
 186     // cells and their amplitude fractions
 187     UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
 188     Double32_t ratios[ncells];
 189
 190     for (Int_t c = 0; c < ncells; c++) {
 191       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
 192
 193       absIds[ncells_true] = digit->GetId();
 194       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 195
 196       if (ratios[ncells_true] > 0.001) ncells_true++;
 197     }
 198
 199     if (ncells_true < 1) {
 200       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
 201       continue;
 202     }
 203
 204     TVector3 gpos;
 205     Float_t g[3];
 206
 207     // calculate new cluster position
 208     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
 209     recPoint->GetGlobalPosition(gpos);
 210     gpos.GetXYZ(g);
 211
 212     // create a new cluster
 213     AliAODCaloCluster *clus = new AliAODCaloCluster();
 214     clus->SetType(AliVCluster::kEMCALClusterv1);
 215     clus->SetE(recPoint->GetEnergy());
 216     clus->SetPosition(g);
 217     clus->SetNCells(ncells_true);
 218     clus->SetCellsAbsId(absIds);
 219     clus->SetCellsAmplitudeFraction(ratios);
 220     // TODO: time not stored
 221     // TODO: some other properties not stored
 222
 223     clusArray->Add(clus);
 224   } // recPoints loop
 225
 226 }
 227
 228 //------------------------------------------------------------------------
 229 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
 230 {
 231   // Unfolds clusters.
 232   //
 233   // Input: TObjArray of clusters, EMCAL cells.
 234   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
 235
 236   Int_t ndigits = 0;
 237
 238   Int_t nclus = clusArray->GetEntriesFast();
 239
 240   /* Fill recPoints with digits
 241   */
 242   for (Int_t i = 0; i < nclus; i++)
 243   {
 244     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
 245     if (!clus->IsEMCAL()) continue;
 246
 247     // new recPoint
 248     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
 249     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 250     fRecPoints->Add(recPoint);
 251
 252     // fill digits
 253     for (Int_t c = 0; c < clus->GetNCells(); c++) {
 254       Int_t absId = clus->GetCellAbsId(c);
 255       Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
 256       Double_t time = cellsEMCAL->GetCellTime(absId);
 257
 258       // NOTE: it is easy to include cells recalibration here:
 259       // amp *= factor;
 260
 261       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
 262
 263       digit->SetId(absId);
 264       digit->SetAmplitude(amp);
 265       digit->SetTime(time);
 266       digit->SetTimeR(time);
 267       digit->SetIndexInList(ndigits);
 268
 269       recPoint->AddDigit(*digit, amp, kFALSE);
 270
 271       ndigits++;
 272     }
 273
 274     // this cluster will be substituted with the result of unfolding
 275     clusArray->RemoveAt(i);
 276     delete clus;
 277   } // cluster loop
 278
 279
 280   /* Peform unfolding
 281   */
 282   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
 283   fClusterUnfolding->MakeUnfolding();
 284
 285   /* Restore clusters from recPoints.
 286   */
 287   RecPoints2Clusters(clusArray);
 288
 289   // clean up
 290   fRecPoints->Delete(); // do not Clear(), it leaks, why?
 291   fDigitsArr->Clear("C");
 292
 293   clusArray->Compress();
 294
 295 }