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