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