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