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