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