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   fRecPoints(NULL),
  84   fDigitsArr(NULL),
  85   fClusterUnfolding(NULL)
  86 {
  87   // Use this constructor, if unsure
  88
  89   Init();
  90 }
  91
  92 //------------------------------------------------------------------------
  93 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
  94   fGeom(NULL),
  95   fLogWeight(logWeight),
  96   fECALocMaxCut(ECALocMaxCut),
  97   fRecPoints(NULL),
  98   fDigitsArr(NULL),
  99   fClusterUnfolding(NULL)
 100 {
 101   // This constructor allows to set parameters
 102   // Recommended values:
 103   //   Float_t logWeight = 4.5, ECALocMaxCut = 0.03
 104
 105   Init();
 106 }
 107
 108 //------------------------------------------------------------------------
 109 void AliEMCALAfterBurnerUF::Init()
 110 {
 111   // After-burner initialization
 112   // Imports geometry.root (if required), creates unfolding class instance
 113   //
 114   // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
 115   //       for this to work, the OCDB geometry must be imported instead
 116
 117   if (!gGeoManager)
 118     Warning("AliEMCALAfterBurnerUF::Init","GeoManager not found, please import the geometry.root file or pass to the geometry the misalignment matrices");
 119 //    TGeoManager::Import("geometry.root");
 120
 121   // required for global cluster position recalculation
 122   if (!gGeoManager)
 123     Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
 124
 125   // initialize geometry, if not yet initialized
 126   if (!AliEMCALGeometry::GetInstance()) {
 127     Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
 128     AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 129   }
 130
 131   // AliEMCALRecPoint is using exactly this call
 132   fGeom = AliEMCALGeometry::GetInstance();
 133   if (!fGeom)
 134     Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
 135
 136   fClusterUnfolding = new AliEMCALUnfolding(fGeom);
 137   fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
 138
 139   // clusters --> recPoints, cells --> digits and back ;)
 140   fRecPoints = new TObjArray(100);
 141   fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
 142 }
 143
 144 //------------------------------------------------------------------------
 145 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
 146 {
 147   if (fClusterUnfolding) delete fClusterUnfolding;
 148
 149   if (fRecPoints) {
 150     fRecPoints->Delete();
 151     delete fRecPoints;
 152   }
 153   if (fDigitsArr) {
 154     fDigitsArr->Clear("C");
 155     delete fDigitsArr;
 156   }
 157 }
 158
 159 //------------------------------------------------------------------------
 160 void AliEMCALAfterBurnerUF::Clear()
 161 {
 162   //Clean the arrays
 163
 164   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
 165   if (fDigitsArr) fDigitsArr->Clear("C");
 166
 167 }
 168 //------------------------------------------------------------------------
 169 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
 170 {
 171   // Restore clusters from recPoints
 172   // Cluster energy, global position, cells and their amplitude fractions are restored
 173   //
 174   // TODO: restore time and other parameters
 175
 176   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
 177   {
 178     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
 179
 180     const Int_t ncells = recPoint->GetMultiplicity();
 181     Int_t ncells_true = 0;
 182
 183     // cells and their amplitude fractions
 184     UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
 185     Double32_t ratios[ncells];
 186
 187     for (Int_t c = 0; c < ncells; c++) {
 188       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
 189
 190       absIds[ncells_true] = digit->GetId();
 191       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 192
 193       if (ratios[ncells_true] > 0.001) ncells_true++;
 194     }
 195
 196     if (ncells_true < 1) {
 197       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
 198       continue;
 199     }
 200
 201     TVector3 gpos;
 202     Float_t g[3];
 203
 204     // calculate new cluster position
 205     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
 206     recPoint->GetGlobalPosition(gpos);
 207     gpos.GetXYZ(g);
 208
 209     // create a new cluster
 210     AliAODCaloCluster *clus = new AliAODCaloCluster();
 211     clus->SetType(AliVCluster::kEMCALClusterv1);
 212     clus->SetE(recPoint->GetEnergy());
 213     clus->SetPosition(g);
 214     clus->SetNCells(ncells_true);
 215     clus->SetCellsAbsId(absIds);
 216     clus->SetCellsAmplitudeFraction(ratios);
 217     // TODO: time not stored
 218     // TODO: some other properties not stored
 219
 220     clusArray->Add(clus);
 221   } // recPoints loop
 222
 223 }
 224
 225 //------------------------------------------------------------------------
 226 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
 227 {
 228   // Unfolds clusters.
 229   //
 230   // Input: TObjArray of clusters, EMCAL cells.
 231   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
 232
 233   Int_t ndigits = 0;
 234
 235   Int_t nclus = clusArray->GetEntriesFast();
 236
 237   /* Fill recPoints with digits
 238   */
 239   for (Int_t i = 0; i < nclus; i++)
 240   {
 241     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
 242     if (!clus->IsEMCAL()) continue;
 243
 244     // new recPoint
 245     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
 246     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 247     fRecPoints->Add(recPoint);
 248
 249     // fill digits
 250     for (Int_t c = 0; c < clus->GetNCells(); c++) {
 251       Int_t absId = clus->GetCellAbsId(c);
 252       Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
 253       Double_t time = cellsEMCAL->GetCellTime(absId);
 254
 255       // NOTE: it is easy to include cells recalibration here:
 256       // amp *= factor;
 257
 258       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
 259
 260       digit->SetId(absId);
 261       digit->SetAmplitude(amp);
 262       digit->SetTime(time);
 263       digit->SetTimeR(time);
 264       digit->SetIndexInList(ndigits);
 265
 266       recPoint->AddDigit(*digit, amp, kFALSE);
 267
 268       ndigits++;
 269     }
 270
 271     // this cluster will be substituted with the result of unfolding
 272     clusArray->RemoveAt(i);
 273     delete clus;
 274   } // cluster loop
 275
 276
 277   /* Peform unfolding
 278   */
 279   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
 280   fClusterUnfolding->MakeUnfolding();
 281
 282   /* Restore clusters from recPoints.
 283   */
 284   RecPoints2Clusters(clusArray);
 285
 286   // clean up
 287   fRecPoints->Delete(); // do not Clear(), it leaks, why?
 288   fDigitsArr->Clear("C");
 289
 290   clusArray->Compress();
 291
 292 }