Adding afterburner class for unfolding studies at the analysis level (Olga)
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Nov 2010 12:24:50 +0000 (12:24 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Nov 2010 12:24:50 +0000 (12:24 +0000)
EMCAL/AliEMCALAfterBurnerUF.cxx [new file with mode: 0644]
EMCAL/AliEMCALAfterBurnerUF.h [new file with mode: 0644]
EMCAL/EMCALrecLinkDef.h
EMCAL/libEMCALrec.pkg

diff --git a/EMCAL/AliEMCALAfterBurnerUF.cxx b/EMCAL/AliEMCALAfterBurnerUF.cxx
new file mode 100644 (file)
index 0000000..5f32d22
--- /dev/null
@@ -0,0 +1,276 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+// After-burner for the EMCAL cluster unfolding algorithm
+//
+// Input: TObjArray  *clusArray -- array of AliAODCaloClusters;
+//        AliAODCaloCells  *cellsEMCAL -- EMCAL cells.
+//
+// Output is appended to clusArray, the original (unfolded or not) clusters
+// are deleted or moved to another position in clusArray.
+//
+// If you want to use particular geometry, you must initialize it _before_
+// creating AliEMCALAfterBurnerUF instance. Add this or similar line to the
+// initialization section:
+//
+//    AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
+//
+// gGeoManager must be initialized for this code to work! Do it yourself or
+// provide geometry.root file in the current directory so that
+// AliEMCALAfterBurnerUF will take it by itself.
+// How to use:
+//
+//   // add this lines to the initialization section of your analysis
+//   AliEMCALAfterBurnerUF *abuf = new AliEMCALAfterBurnerUF();
+//   TObjArray *clusArray = new TObjArray(100);
+//
+//
+//   AliAODEvent *event = (AliAODEvent*) InputEvent();
+//   AliAODCaloCells *cellsEMCAL = event->GetEMCALCells();
+//
+//   for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
+//   {
+//     AliAODCaloCluster *clus = event->GetCaloCluster(i);
+//
+//     clusArray->Add(clus->Clone());   // NOTE _CLONE_ in this line
+//   }
+//
+//   abuf->UnfoldClusters(clusArray, cellsEMCAL);
+//
+//   // do an analysis with clusArray
+//   // ....
+//
+//   // prevent memory leak
+//   clusArray->Delete();
+//
+//----
+//  Author: Olga Driga (SUBATECH)
+
+// --- ROOT system ---
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TGeoManager.h>
+
+// --- Standard library --
+
+// --- AliRoot header files ---
+#include <AliEMCALAfterBurnerUF.h>
+#include <AliEMCALGeometry.h>
+#include <AliEMCALUnfolding.h>
+#include <AliAODCaloCluster.h>
+#include <AliAODCaloCells.h>
+#include <AliEMCALRecPoint.h>
+#include <AliEMCALDigit.h>
+
+ClassImp(AliEMCALAfterBurnerUF)
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
+  fGeom(NULL),
+  fLogWeight(4.5),      // correct?
+  fECALocMaxCut(0.03),  // value suggested by Adam
+  fRecPoints(NULL),
+  fDigitsArr(NULL),
+  fClusterUnfolding(NULL)
+{
+  // Use this constructor, if unsure
+
+  Init();
+}
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
+  fGeom(NULL),
+  fLogWeight(logWeight),
+  fECALocMaxCut(ECALocMaxCut),
+  fRecPoints(NULL),
+  fDigitsArr(NULL),
+  fClusterUnfolding(NULL)
+{
+  // This constructor allows to set parameters
+  // Recommended values:
+  //   Float_t logWeight = 4.5, ECALocMaxCut = 0.03
+
+  Init();
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::Init()
+{
+  // After-burner initialization
+  // Imports geometry.root (if required), creates unfolding class instance
+  //
+  // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
+  //       for this to work, the OCDB geometry must be imported instead
+
+  if (!gGeoManager)
+    TGeoManager::Import("geometry.root");
+
+  // required for global cluster position recalculation
+  if (!gGeoManager)
+    Fatal("AliEMCALAfterBurnerUF::Init", "failed to import geometry.root");
+
+  // initialize geometry, if not yet initialized
+  if (!AliEMCALGeometry::GetInstance()) {
+    Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETE");
+    AliEMCALGeometry::GetInstance("EMCAL_COMPLETE");
+  }
+
+  // AliEMCALRecPoint is using exactly this call
+  fGeom = AliEMCALGeometry::GetInstance();
+  if (!fGeom)
+    Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
+
+  fClusterUnfolding = new AliEMCALUnfolding(fGeom);
+  fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
+
+  // clusters --> recPoints, cells --> digits and back ;)
+  fRecPoints = new TObjArray(100);
+  fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
+}
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
+{
+  if (fClusterUnfolding) delete fClusterUnfolding;
+
+  if (fRecPoints) delete fRecPoints;
+  if (fDigitsArr) delete fDigitsArr;
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
+{
+  // Restore clusters from recPoints
+  // Cluster energy, global position, cells and their amplitude fractions are restored
+  //
+  // TODO: restore time and other parameters
+
+  for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
+  {
+    AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
+
+    Int_t ncells = recPoint->GetMultiplicity();
+    Int_t ncells_true = 0;
+
+    // cells and their amplitude fractions
+    UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
+    Double32_t ratios[ncells];
+
+    for (Int_t c = 0; c < ncells; c++) {
+      AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
+
+      absIds[ncells_true] = digit->GetId();
+      ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
+
+      if (ratios[ncells_true] > 0.001) ncells_true++;
+    }
+
+    if (ncells_true < 1) {
+      Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
+      continue;
+    }
+
+    TVector3 gpos;
+    Float_t g[3];
+
+    // calculate new cluster position
+    recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
+    recPoint->GetGlobalPosition(gpos);
+    gpos.GetXYZ(g);
+
+    // create a new cluster
+    AliAODCaloCluster *clus = new AliAODCaloCluster();
+    clus->SetType(AliAODCaloCluster::kEMCALClusterv1);
+    clus->SetE(recPoint->GetEnergy());
+    clus->SetPosition(g);
+    clus->SetNCells(ncells_true);
+    clus->SetCellsAbsId(absIds);
+    clus->SetCellsAmplitudeFraction(ratios);
+    // TODO: time not stored
+    // TODO: some other properties not stored
+
+    clusArray->Add(clus);
+  } // recPoints loop
+
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliAODCaloCells *cellsEMCAL)
+{
+  // Unfolds clusters.
+  //
+  // Input: TObjArray of clusters, EMCAL cells.
+  // Output is added to the same array, original clusters are _deleted_ or moved to another position.
+
+  Int_t ndigits = 0;
+
+  Int_t nclus = clusArray->GetEntriesFast();
+
+  /* Fill recPoints with digits
+  */
+  for (Int_t i = 0; i < nclus; i++)
+  {
+    AliAODCaloCluster *clus = (AliAODCaloCluster*) clusArray->At(i);
+    if (!clus->IsEMCAL()) continue;
+
+    // new recPoint
+    AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
+    recPoint->SetClusterType(AliAODCaloCluster::kEMCALClusterv1);
+    fRecPoints->Add(recPoint);
+
+    // fill digits
+    for (Int_t c = 0; c < clus->GetNCells(); c++) {
+      Int_t absId = clus->GetCellAbsId(c);
+      Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
+      Double_t time = cellsEMCAL->GetCellTime(absId);
+
+      // NOTE: it is easy to include cells recalibration here:
+      // amp *= factor;
+
+      AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
+
+      digit->SetId(absId);
+      digit->SetAmplitude(amp);
+      digit->SetTime(time);
+      digit->SetTimeR(time);
+      digit->SetIndexInList(ndigits);
+
+      recPoint->AddDigit(*digit, amp, kFALSE);
+
+      ndigits++;
+    }
+
+    // this cluster will be substituted with the result of unfolding
+    clusArray->RemoveAt(i);
+    delete clus;
+  } // cluster loop
+
+
+  /* Peform unfolding
+  */
+  fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
+  fClusterUnfolding->MakeUnfolding();
+
+  /* Restore clusters from recPoints.
+  */
+  RecPoints2Clusters(clusArray);
+
+  // clean up
+  fRecPoints->Delete();
+  fDigitsArr->Delete();
+
+  clusArray->Compress();
+
+}
diff --git a/EMCAL/AliEMCALAfterBurnerUF.h b/EMCAL/AliEMCALAfterBurnerUF.h
new file mode 100644 (file)
index 0000000..b3c81ce
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALIEMCALAFTERBURNERUF_H
+#define ALIEMCALAFTERBURNERUF_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//_________________________________________________________________________
+//  After-burner for the EMCAL cluster unfolding algorithm
+//
+//  See cxx for details on how to use it
+//
+//  Author: Olga Driga (SUBATECH)
+//
+
+// --- ROOT system ---
+class TObjArray;
+class TClonesArray;
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+class AliEMCALGeometry;
+class AliEMCALUnfolding;
+class AliAODCaloCells;
+
+class AliEMCALAfterBurnerUF {
+
+  public:
+    AliEMCALAfterBurnerUF();
+    AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut);
+    virtual ~AliEMCALAfterBurnerUF();
+
+  private:
+    AliEMCALAfterBurnerUF(const AliEMCALAfterBurnerUF & uf) ; // cpy ctor not needed, put here to avoid compilation warning 
+    AliEMCALAfterBurnerUF & operator = (const AliEMCALAfterBurnerUF & uf) ;//cpy assignment, put here to avoid compilation warning 
+  
+  public:
+    virtual void Init();
+    virtual void RecPoints2Clusters(TObjArray *clusArray);
+    virtual void UnfoldClusters(TObjArray *clusArray, AliAODCaloCells *cellsEMCAL);  // does the job
+
+    // getters and setters
+    virtual AliEMCALUnfolding *GetClusterUnfoldingInstance() { return fClusterUnfolding; }
+
+  protected:
+    AliEMCALGeometry  *fGeom;          // EMCAL geometry
+    Float_t            fLogWeight;     // used in AliEMCALRecPoint::EvalGlobalPosition()
+    Float_t            fECALocMaxCut;  // this amount of energy must distinguish a local maximum from its neighbours
+    TObjArray         *fRecPoints;     // cluster <=> recPoint
+    TClonesArray      *fDigitsArr;     // cell <=> digit
+
+    AliEMCALUnfolding *fClusterUnfolding;  // unfolding class instance
+
+    ClassDef(AliEMCALAfterBurnerUF,1)
+} ;
+
+#endif // AliEMCALAFTERBURNERUF_H
index 5ca5b2c..f010c07 100644 (file)
@@ -12,7 +12,7 @@
 #pragma link C++ class AliEMCALPID+;
 #pragma link C++ class AliEMCALQADataMakerRec+;
 #pragma link C++ class AliEMCALAodCluster+;
-
+#pragma link C++ class AliEMCALAfterBurnerUF+;
 
 
 #endif
index d4d9c10..4eec35a 100644 (file)
@@ -10,8 +10,8 @@ AliEMCALUnfolding.cxx \
 AliEMCALTracker.cxx \
 AliEMCALPID.cxx \
 AliEMCALQADataMakerRec.cxx \
-AliEMCALAodCluster.cxx
-
+AliEMCALAodCluster.cxx \
+AliEMCALAfterBurnerUF.cxx
 
 HDRS= $(SRCS:.cxx=.h)