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   fIsInputCalibrated(kFALSE),
  65   fJustClusters(kFALSE),
  66   fDigitsArr(NULL),
  67   fTreeR(NULL),
  68   fRecPoints(NULL),
  69   fGeom(NULL),
  70   fCalibData(NULL),
  71   fCaloPed(NULL),
  72   fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
  73   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  74   fDefaultInit(kFALSE),fToUnfold(kFALSE),
  75   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  76   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
  77   fClusterUnfolding(NULL)
  78 {
  79   // ctor
  80
  81   Init();
  82 }
  83
  84 //____________________________________________________________________________
  85 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
  86   fIsInputCalibrated(kFALSE),
  87   fJustClusters(kFALSE),
  88   fDigitsArr(NULL),
  89   fTreeR(NULL),
  90   fRecPoints(NULL),
  91   fGeom(geometry),
  92   fCalibData(NULL),
  93   fCaloPed(NULL),
  94   fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
  95   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  96   fDefaultInit(kFALSE),fToUnfold(kFALSE),
  97   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  98   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
  99   fClusterUnfolding(NULL)
 100 {
 101   // Ctor with the indication of the file where header Tree and digits Tree are stored.
 102   // Use this contructor to avoid usage of Init() which uses runloader.
 103   // Change needed by HLT - MP.
 104   // Note for the future: the use on runloader should be avoided or optional at least.
 105   // Another way is to make Init virtual and protected at least
 106   // such that the deriving classes can overload Init();
 107
 108   if (!fGeom)
 109   {
 110     AliFatal("Geometry not initialized.");
 111   }
 112   Int_t i=0;
 113   for (i = 0; i < 8; i++)
 114     fSSPars[i] = 0.;
 115   for (i = 0; i < 3; i++) {
 116     fPar5[i] = 0.;
 117     fPar6[i] = 0.;
 118   }
 119 }
 120
 121 //____________________________________________________________________________
 122 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry,
 123                                          AliEMCALCalibData *calib,
 124                                          AliCaloCalibPedestal *caloped):
 125   fIsInputCalibrated(kFALSE),
 126   fJustClusters(kFALSE),
 127   fDigitsArr(NULL),
 128   fTreeR(NULL),
 129   fRecPoints(NULL),
 130   fGeom(geometry),
 131   fCalibData(calib),
 132   fCaloPed(caloped),
 133   fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
 134   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
 135   fDefaultInit(kFALSE),fToUnfold(kFALSE),
 136   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
 137   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
 138   fClusterUnfolding(NULL)
 139 {
 140   // ctor, geometry and calibration are initialized elsewhere.
 141
 142   if (!fGeom)
 143     AliFatal("Geometry not initialized.");
 144
 145   Int_t i=0;
 146   for (i = 0; i < 8; i++)
 147     fSSPars[i] = 0.;
 148   for (i = 0; i < 3; i++) {
 149     fPar5[i] = 0.;
 150     fPar6[i] = 0.;
 151   }
 152 }
 153
 154 //____________________________________________________________________________
 155 AliEMCALClusterizer::~AliEMCALClusterizer()
 156 {
 157   // dtor
 158   //Already deleted in AliEMCALReconstructor.
 159
 160   if(fClusterUnfolding) delete fClusterUnfolding;
 161
 162   // make sure we delete the rec points array
 163   DeleteRecPoints();
 164
 165   //Delete digits array
 166   DeleteDigits();
 167
 168 }
 169
 170 //____________________________________________________________________________
 171 void AliEMCALClusterizer::DeleteRecPoints()
 172 {
 173   // free the cluster array
 174   if (fRecPoints)
 175     {
 176       AliDebug(2, "Deleting fRecPoints.");
 177       fRecPoints->Delete();
 178       delete fRecPoints;
 179       fRecPoints = 0;
 180     }
 181 }
 182
 183 //____________________________________________________________________________
 184 void AliEMCALClusterizer::DeleteDigits()
 185 {
 186   // free the digits array
 187   if (fDigitsArr)
 188     {
 189       AliDebug(2, "Deleting fDigitsArr.");
 190       fDigitsArr->Clear("C");
 191       delete fDigitsArr;
 192       fDigitsArr = 0;
 193     }
 194 }
 195
 196 //____________________________________________________________________________
 197 void AliEMCALClusterizer::Calibrate(Float_t & amp, Float_t & time, const Int_t absId)
 198 {
 199   // Convert digitized amplitude into energy, calibrate time
 200   // Calibration parameters are taken from OCDB : OCDB/EMCAL/Calib/Data
 201
 202   //Check if time is too large or too small, indication of a noisy channel, remove in this case
 203   if(time > fTimeMax || time < fTimeMin) {
 204     amp  = 0 ;
 205     time = 0 ;
 206     return ;
 207   }
 208
 209   //Return energy with default parameters if calibration is not available
 210   if (!fCalibData && !fCaloPed) {
 211     if (fIsInputCalibrated == kTRUE)
 212     {
 213       AliDebug(10, Form("Input already calibrated!"));
 214       return ;
 215     }
 216     else{
 217       AliFatal("OCDB calibration and bad map parameters are not available");
 218       return;
 219     }
 220   }
 221
 222   if (fGeom==0)
 223     AliFatal("Did not get geometry from EMCALLoader") ;
 224
 225   Int_t iSupMod = -1;
 226   Int_t nModule = -1;
 227   Int_t nIphi   = -1;
 228   Int_t nIeta   = -1;
 229   Int_t iphi    = -1;
 230   Int_t ieta    = -1;
 231
 232   Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 233   if(!bCell) {
 234     fGeom->PrintGeometry();
 235     AliError(Form("Wrong cell id number : %i", absId));
 236     //assert(0); // GCB: This aborts reconstruction of raw simulations
 237     //where simulation had more SM than default geometry,
 238     //change to return 0, to avoid aborting good generations.
 239     amp  = 0;
 240     time = 0;
 241     return ;
 242   }
 243
 244   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 245
 246   // Check if channel is bad (dead or hot), in this case return 0.
 247   // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
 248   // for the moment keep it here but remember to do the selection at the sdigitizer level
 249   // and remove it from here
 250   if (fCaloPed) {
 251     Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
 252     if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
 253       AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
 254       amp  = 0 ;
 255       time = 0 ;
 256       return ;
 257     }
 258   }
 259
 260   if (fIsInputCalibrated || !fCalibData)
 261   {
 262     AliDebug(10, Form("Input already calibrated!"));
 263     return ;
 264   }
 265
 266   fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 267   fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 268   fTimeECA        = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
 269
 270   time -= fTimeECA ;
 271   amp   = amp * fADCchannelECA - fADCpedestalECA ;
 272
 273 }
 274
 275 //____________________________________________________________________________
 276 void AliEMCALClusterizer::GetCalibrationParameters()
 277 {
 278   // Set calibration parameters:
 279   // If calibration database exists, they are read from database,
 280   // otherwise, they are taken from digitizer.
 281   // It is a user responsilibity to open CDB before reconstruction,
 282   // for example:
 283   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 284
 285   if (fIsInputCalibrated)
 286     return;
 287
 288   //Check if calibration is stored in data base
 289   if(!fCalibData)
 290   {
 291     AliCDBEntry *entry = (AliCDBEntry*)
 292     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 293     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 294   }
 295
 296   if(!fCalibData)
 297     AliFatal("Calibration parameters not found in CDB!");
 298 }
 299
 300 //____________________________________________________________________________
 301 void AliEMCALClusterizer::GetCaloCalibPedestal()
 302 {
 303   // Set calibration parameters:
 304   // if calibration database exists, they are read from database,
 305   // otherwise, they are taken from digitizer.
 306   //
 307   // It is a user responsilibity to open CDB before reconstruction,
 308   // for example:
 309   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 310
 311   if (fIsInputCalibrated)
 312     return;
 313
 314   // Check if calibration is stored in data base
 315   if(!fCaloPed)
 316     {
 317       AliCDBEntry *entry = (AliCDBEntry*)
 318         AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 319       if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
 320     }
 321
 322   if(!fCaloPed)
 323     AliFatal("Pedestal info not found in CDB!");
 324 }
 325
 326 //____________________________________________________________________________
 327 void AliEMCALClusterizer::Init()
 328 {
 329   // Make all memory allocations which can not be done in default constructor.
 330   // Attach the Clusterizer task to the list of EMCAL tasks
 331
 332   AliRunLoader *rl = AliRunLoader::Instance();
 333   if (rl->GetAliRun()){
 334     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
 335     if(emcal)fGeom = emcal->GetGeometry();
 336   }
 337
 338   if(!fGeom){
 339     fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
 340   }
 341
 342   AliDebug(1,Form("geom %p",fGeom));
 343
 344   if(!gMinuit)
 345     gMinuit = new TMinuit(100) ;
 346
 347   Int_t i=0;
 348   for (i = 0; i < 8; i++)
 349     fSSPars[i] = 0.;
 350   for (i = 0; i < 3; i++) {
 351     fPar5[i] = 0.;
 352     fPar6[i] = 0.;
 353   }
 354 }
 355
 356 //____________________________________________________________________________
 357 void AliEMCALClusterizer::InitParameters()
 358 {
 359   // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
 360
 361   return InitParameters(AliEMCALReconstructor::GetRecParam());
 362 }
 363
 364 //____________________________________________________________________________
 365 void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
 366 {
 367   // Initializes the parameters for the Clusterizer
 368
 369   fNumberOfECAClusters = 0 ;
 370   fCalibData           = 0 ;
 371   fCaloPed             = 0 ;
 372
 373   if(!recParam) {
 374     AliFatal("Reconstruction parameters for EMCAL not set!");
 375   }
 376
 377   fECAClusteringThreshold = recParam->GetClusteringThreshold();
 378   fECAW0                  = recParam->GetW0();
 379   fMinECut                = recParam->GetMinECut();
 380   fToUnfold               = recParam->GetUnfold();
 381   fECALocMaxCut           = recParam->GetLocMaxCut();
 382   fTimeCut                = recParam->GetTimeCut();
 383   fTimeMin                = recParam->GetTimeMin();
 384   fTimeMax                = recParam->GetTimeMax();
 385
 386   //For NxN
 387   SetNRowDiff(recParam->GetNRowDiff());
 388   SetNColDiff(recParam->GetNColDiff());
 389
 390   AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
 391                   "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
 392                   fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
 393
 394   if (fToUnfold) {
 395     Int_t i=0;
 396     for (i = 0; i < 8; i++) {
 397       fSSPars[i] = recParam->GetSSPars(i);
 398     } //end of loop over parameters
 399     for (i = 0; i < 3; i++) {
 400       fPar5[i] = recParam->GetPar5(i);
 401       fPar6[i] = recParam->GetPar6(i);
 402     } //end of loop over parameters
 403
 404     InitClusterUnfolding();
 405
 406     for (i = 0; i < 8; i++) {
 407       AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
 408     }
 409     for (i = 0; i < 3; i++) {
 410       AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
 411       AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
 412     }
 413   } // to unfold
 414 }
 415
 416 //____________________________________________________________________________
 417 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
 418 {
 419   // Print clusterizer parameters
 420
 421   TString message("\n") ;
 422
 423   if (strcmp(GetName(),"") == 0) {
 424     printf("AliEMCALClusterizer not initialized\n");
 425     return;
 426   }
 427
 428   // Print parameters
 429   TString taskName(Version()) ;
 430
 431   printf("--------------- ");
 432   printf("%s",taskName.Data()) ;
 433   printf(" ");
 434   printf("Clusterizing digits: ");
 435   printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 436   printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 437   if (fToUnfold) {
 438     printf("\nUnfolding on\n");
 439     printf("Unfolding parameters: fSSpars: \n");
 440     Int_t i=0;
 441     for (i = 0; i < 8; i++) {
 442       printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
 443     }
 444     printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
 445     for (i = 0; i < 3; i++) {
 446       printf("fPar5[%d] = %f \n", i, fPar5[i]);
 447       printf("fPar6[%d] = %f \n", i, fPar6[i]);
 448     }
 449   }
 450   else
 451     printf("\nUnfolding off\n");
 452
 453   printf("------------------------------------------------------------------");
 454 }
 455
 456 //____________________________________________________________________________
 457 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
 458 {
 459   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 460
 461   if (strstr(option,"deb")) {
 462     printf("PrintRecPoints: Clusterization result:") ;
 463     printf("           Found %d ECA Rec Points\n ",
 464            fRecPoints->GetEntriesFast()) ;
 465   }
 466
 467   if (strstr(option,"all")) {
 468     if (strstr(option,"deb")) {
 469       printf("\n-----------------------------------------------------------------------\n") ;
 470       printf("Clusters in ECAL section\n") ;
 471       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 472     }
 473     Int_t index;
 474     for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
 475       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
 476       if (!rp)
 477         continue;
 478
 479       TVector3  globalpos;
 480       //rp->GetGlobalPosition(globalpos);
 481       TVector3  localpos;
 482       rp->GetLocalPosition(localpos);
 483       Float_t lambda[2];
 484       rp->GetElipsAxis(lambda);
 485
 486       Int_t nprimaries=0;
 487       Int_t * primaries = rp->GetPrimaries(nprimaries);
 488       if(strstr(option,"deb"))
 489         printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 490                rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 491                globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 492                rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 493       if(strstr(option,"deb")){
 494         for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 495           printf("%d ", primaries[iprimary] ) ;
 496         }
 497       }
 498     }
 499
 500     if(strstr(option,"deb"))
 501       printf("\n-----------------------------------------------------------------------\n");
 502   }
 503 }
 504
 505 //___________________________________________________________________
 506 void  AliEMCALClusterizer::PrintRecoInfo()
 507 {
 508   // Print reco version
 509
 510   printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
 511 }
 512
 513 //____________________________________________________________________________
 514 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
 515 {
 516   // Read the digits from the input tree
 517
 518   TBranch *branch = digitsTree->GetBranch("EMCAL");
 519   if (!branch) {
 520     AliError("can't get the branch with the EMCAL digits !");
 521     return;
 522   }
 523   if (!fDigitsArr)
 524     fDigitsArr = new TClonesArray("AliEMCALDigit",100);
 525   branch->SetAddress(&fDigitsArr);
 526   branch->GetEntry(0);
 527 }
 528
 529 //____________________________________________________________________________
 530 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
 531 {
 532   // Read the digits from the input tree
 533
 534   AliDebug(9, "Making array for EMCAL clusters");
 535   fRecPoints = new TObjArray(1000);
 536   if (clustersTree) {
 537     fTreeR = clustersTree;
 538     Int_t split = 0;
 539     Int_t bufsize = 32000;
 540     fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
 541   }
 542 }
 543
 544 //___________________________________________________________________
 545 void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
 546 {
 547   // Flag to indicate that input is calibrated - the case when we run already on ESD
 548
 549   fIsInputCalibrated = val;
 550 }
 551
 552 //___________________________________________________________________
 553 void AliEMCALClusterizer::SetJustClusters(Bool_t val)
 554 {
 555   // Flag to indicate that we are running on ESDs, when calling
 556   // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes
 557
 558   fJustClusters = val;
 559 }