EMCAL/AliEMCALClusterizer.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
  16 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 //  Base class for the clusterization algorithm (pure abstract)
  20 //*--
  21 //*-- Author: Yves Schutz  SUBATECH
  22 //
  23 //   Clusterization mother class. Contains common methods/data members of different
  24 //   clusterizers. GCB 2010
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27 // --- ROOT system ---
  28 #include "TClonesArray.h"
  29 #include "TTree.h"
  30 #include <TFile.h>
  31 class TFolder;
  32 #include <TMath.h>
  33 #include <TMinuit.h>
  34 #include <TTree.h>
  35 class TSystem;
  36 #include <TBenchmark.h>
  37 #include <TBrowser.h>
  38 #include <TROOT.h>
  39
  40 // --- Standard library ---
  41 #include <cassert>
  42
  43 // --- AliRoot header files ---
  44 #include "AliEMCALClusterizer.h"
  45 #include "AliEMCALReconstructor.h"
  46 #include "AliRunLoader.h"
  47 #include "AliRun.h"
  48 #include "AliLog.h"
  49 #include "AliEMCAL.h"
  50 #include "AliEMCALRecPoint.h"
  51 #include "AliEMCALRecParam.h"
  52 #include "AliEMCALGeometry.h"
  53 #include "AliEMCALRecParam.h"
  54 #include "AliCDBManager.h"
  55 #include "AliCaloCalibPedestal.h"
  56 #include "AliEMCALCalibData.h"
  57 class AliCDBStorage;
  58 #include "AliCDBEntry.h"
  59
  60 ClassImp(AliEMCALClusterizer)
  61
  62 //____________________________________________________________________________
  63 AliEMCALClusterizer::AliEMCALClusterizer():
  64   fDigitsArr(NULL),
  65   fTreeR(NULL),
  66   fRecPoints(NULL),
  67   fGeom(NULL),
  68   fCalibData(NULL),
  69   fCaloPed(NULL),
  70   fADCchannelECA(0.),fADCpedestalECA(0.),
  71   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  72   fDefaultInit(kFALSE),fToUnfold(kFALSE),
  73   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  74   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
  75 {
  76   // ctor
  77
  78   Init();
  79
  80 }
  81
  82 //____________________________________________________________________________
  83 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
  84   fDigitsArr(NULL),
  85   fTreeR(NULL),
  86   fRecPoints(NULL),
  87   fGeom(geometry),
  88   fCalibData(NULL),
  89   fCaloPed(NULL),
  90   fADCchannelECA(0.),fADCpedestalECA(0.),
  91   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  92   fDefaultInit(kFALSE),fToUnfold(kFALSE),
  93   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  94   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
  95 {
  96   // ctor with the indication of the file where header Tree and digits Tree are stored
  97   // use this contructor to avoid usage of Init() which uses runloader
  98   // change needed by HLT - MP
  99
 100   // Note for the future: the use on runloader should be avoided or optional at least
 101   // another way is to make Init virtual and protected at least such that the deriving classes can overload
 102   // Init() ;
 103   //
 104
 105   if (!fGeom)
 106   {
 107     AliFatal("Geometry not initialized.");
 108   }
 109
 110 }
 111
 112 //____________________________________________________________________________
 113 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped):
 114   fDigitsArr(NULL),
 115   fTreeR(NULL),
 116   fRecPoints(NULL),
 117   fGeom(geometry),
 118   fCalibData(calib),
 119   fCaloPed(caloped),
 120   fADCchannelECA(0.),fADCpedestalECA(0.),
 121   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
 122   fDefaultInit(kFALSE),fToUnfold(kFALSE),
 123   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
 124   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
 125 {
 126         // ctor, geometry and calibration are initialized elsewhere.
 127
 128         if (!fGeom)
 129                 AliFatal("Geometry not initialized.");
 130
 131 }
 132
 133
 134 //____________________________________________________________________________
 135 AliEMCALClusterizer::~AliEMCALClusterizer()
 136 {
 137   // dtor
 138   //Already deleted in AliEMCALReconstructor.
 139
 140 //   if(fGeom)      delete fGeom;
 141 //   if(fCalibData) delete fCalibData;
 142 //   if(fCaloPed)   delete fCaloPed;
 143
 144 //   if (fDigitsArr) {
 145 //     fDigitsArr->Clear("C");
 146 //     delete fDigitsArr;
 147 //   }
 148 //   if (fRecPoints) {
 149 //     fRecPoints->Delete();
 150 //     delete fRecPoints;
 151 //  }
 152 }
 153
 154 //____________________________________________________________________________
 155 Float_t  AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId)
 156 {
 157
 158   // Convert digitized amplitude into energy.
 159   // Calibration parameters are taken from calibration data base for raw data,
 160   // or from digitizer parameters for simulated data.
 161   if(fCalibData){
 162
 163     if (fGeom==0)
 164       AliFatal("Did not get geometry from EMCALLoader") ;
 165
 166     Int_t iSupMod = -1;
 167     Int_t nModule = -1;
 168     Int_t nIphi   = -1;
 169     Int_t nIeta   = -1;
 170     Int_t iphi    = -1;
 171     Int_t ieta    = -1;
 172
 173     Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 174     if(!bCell) {
 175       fGeom->PrintGeometry();
 176       Error("Calibrate()"," Wrong cell id number : %i", absId);
 177       assert(0);
 178     }
 179
 180     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 181
 182     // Check if channel is bad (dead or hot), in this case return 0.
 183     // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
 184     // for the moment keep it here but remember to do the selection at the sdigitizer level
 185     // and remove it from here
 186     Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
 187     if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
 188                   AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
 189                   return 0;
 190     }
 191     //Check if time is too large or too small, indication of a noisy channel, remove in this case
 192     if(time > fTimeMax || time < fTimeMin) return 0;
 193
 194     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 195     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 196
 197     return -fADCpedestalECA + amp * fADCchannelECA ;
 198
 199   }
 200   else //Return energy with default parameters if calibration is not available
 201     return -fADCpedestalECA + amp * fADCchannelECA ;
 202
 203 }
 204
 205 //____________________________________________________________________________
 206 void AliEMCALClusterizer::GetCalibrationParameters()
 207 {
 208   // Set calibration parameters:
 209   // if calibration database exists, they are read from database,
 210   // otherwise, they are taken from digitizer.
 211   //
 212   // It is a user responsilibity to open CDB before reconstruction,
 213   // for example:
 214   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 215
 216   //Check if calibration is stored in data base
 217
 218   if(!fCalibData)
 219   {
 220     AliCDBEntry *entry = (AliCDBEntry*)
 221     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 222     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 223   }
 224
 225   if(!fCalibData)
 226     AliFatal("Calibration parameters not found in CDB!");
 227
 228 }
 229
 230 //____________________________________________________________________________
 231 void AliEMCALClusterizer::GetCaloCalibPedestal()
 232 {
 233         // Set calibration parameters:
 234         // if calibration database exists, they are read from database,
 235         // otherwise, they are taken from digitizer.
 236         //
 237         // It is a user responsilibity to open CDB before reconstruction,
 238         // for example:
 239         // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 240
 241         //Check if calibration is stored in data base
 242
 243         if(!fCaloPed)
 244   {
 245                 AliCDBEntry *entry = (AliCDBEntry*)
 246                 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 247                 if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
 248   }
 249
 250         if(!fCaloPed)
 251                 AliFatal("Pedestal info not found in CDB!");
 252
 253 }
 254
 255 //____________________________________________________________________________
 256 void AliEMCALClusterizer::Init()
 257 {
 258   // Make all memory allocations which can not be done in default constructor.
 259   // Attach the Clusterizer task to the list of EMCAL tasks
 260
 261   AliRunLoader *rl = AliRunLoader::Instance();
 262   if (rl->GetAliRun() && dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL")))
 263     fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
 264   else
 265     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
 266
 267   AliDebug(1,Form("geom %p",fGeom));
 268
 269   if(!gMinuit)
 270     gMinuit = new TMinuit(100) ;
 271
 272 }
 273
 274 //____________________________________________________________________________
 275 void AliEMCALClusterizer::InitParameters()
 276 {
 277   // Initializes the parameters for the Clusterizer
 278   fNumberOfECAClusters = 0 ;
 279   fCalibData           = 0 ;
 280   fCaloPed             = 0 ;
 281
 282   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 283   if(!recParam) {
 284     AliFatal("Reconstruction parameters for EMCAL not set!");
 285   } else {
 286     fECAClusteringThreshold = recParam->GetClusteringThreshold();
 287     fECAW0                  = recParam->GetW0();
 288     fMinECut                = recParam->GetMinECut();
 289     fToUnfold               = recParam->GetUnfold();
 290     if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");
 291     fECALocMaxCut           = recParam->GetLocMaxCut();
 292     fTimeCut                = recParam->GetTimeCut();
 293     fTimeMin                = recParam->GetTimeMin();
 294     fTimeMax                = recParam->GetTimeMax();
 295
 296     AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
 297                     fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
 298   }
 299
 300 }
 301
 302
 303 //____________________________________________________________________________
 304 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
 305 {
 306   // Print clusterizer parameters
 307
 308   TString message("\n") ;
 309
 310   if( strcmp(GetName(), "") !=0 ){
 311
 312     // Print parameters
 313
 314     TString taskName(Version()) ;
 315
 316     printf("--------------- ");
 317     printf("%s",taskName.Data()) ;
 318     printf(" ");
 319     printf("Clusterizing digits: ");
 320     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 321     printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 322     if(fToUnfold)
 323       printf("\nUnfolding on\n");
 324     else
 325       printf("\nUnfolding off\n");
 326
 327     printf("------------------------------------------------------------------");
 328   }
 329   else
 330     printf("AliEMCALClusterizer not initialized ") ;
 331 }
 332
 333 //____________________________________________________________________________
 334 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
 335 {
 336   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 337   if(strstr(option,"deb")) {
 338     printf("PrintRecPoints: Clusterization result:") ;
 339
 340     printf("           Found %d ECA Rec Points\n ",
 341            fRecPoints->GetEntriesFast()) ;
 342   }
 343
 344   if(strstr(option,"all")) {
 345     if(strstr(option,"deb")) {
 346       printf("\n-----------------------------------------------------------------------\n") ;
 347       printf("Clusters in ECAL section\n") ;
 348       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 349     }
 350     Int_t index;
 351     for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
 352       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
 353       if(rp){
 354         TVector3  globalpos;
 355         //rp->GetGlobalPosition(globalpos);
 356         TVector3  localpos;
 357         rp->GetLocalPosition(localpos);
 358         Float_t lambda[2];
 359         rp->GetElipsAxis(lambda);
 360
 361         Int_t nprimaries=0;
 362         Int_t * primaries = rp->GetPrimaries(nprimaries);
 363
 364         if(strstr(option,"deb"))
 365           printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 366                  rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 367                  globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 368                  rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 369         if(strstr(option,"deb")){
 370           for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 371             printf("%d ", primaries[iprimary] ) ;
 372           }
 373         }
 374       }
 375     }
 376
 377     if(strstr(option,"deb"))
 378       printf("\n-----------------------------------------------------------------------\n");
 379   }
 380 }
 381
 382 //___________________________________________________________________
 383 void  AliEMCALClusterizer::PrintRecoInfo()
 384 {
 385   //print reco version
 386   printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
 387
 388 }
 389
 390 //____________________________________________________________________________
 391 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
 392 {
 393   // Read the digits from the input tree
 394   TBranch *branch = digitsTree->GetBranch("EMCAL");
 395   if (!branch) {
 396     AliError("can't get the branch with the EMCAL digits !");
 397     return;
 398   }
 399   if (!fDigitsArr)
 400     fDigitsArr = new TClonesArray("AliEMCALDigit",100);
 401   branch->SetAddress(&fDigitsArr);
 402   branch->GetEntry(0);
 403 }
 404
 405 //____________________________________________________________________________
 406 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
 407 {
 408   // Read the digits from the input tree
 409   fTreeR = clustersTree;
 410
 411   AliDebug(9, "Making array for EMCAL clusters");
 412   fRecPoints = new TObjArray(100) ;
 413   Int_t split = 0;
 414   Int_t bufsize = 32000;
 415   fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
 416 }
 417
 418
 419