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