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 //____________________________________________________________________________
  64 AliEMCALClusterizer::AliEMCALClusterizer():
  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 //____________________________________________________________________________
  85 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
  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
 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 such that the deriving classes can overload
 105   // Init() ;
 106   //
 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, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped):
 123   fDigitsArr(NULL),
 124   fTreeR(NULL),
 125   fRecPoints(NULL),
 126   fGeom(geometry),
 127   fCalibData(calib),
 128   fCaloPed(caloped),
 129   fADCchannelECA(0.),fADCpedestalECA(0.),
 130   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
 131   fDefaultInit(kFALSE),fToUnfold(kFALSE),
 132   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
 133   fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
 134   fClusterUnfolding(NULL)
 135 {
 136   // ctor, geometry and calibration are initialized elsewhere.
 137
 138   if (!fGeom)
 139     AliFatal("Geometry not initialized.");
 140
 141   Int_t i=0;
 142   for (i = 0; i < 8; i++)
 143     fSSPars[i] = 0.;
 144   for (i = 0; i < 3; i++) {
 145     fPar5[i] = 0.;
 146     fPar6[i] = 0.;
 147   }
 148
 149 }
 150
 151
 152 //____________________________________________________________________________
 153 AliEMCALClusterizer::~AliEMCALClusterizer()
 154 {
 155   // dtor
 156   //Already deleted in AliEMCALReconstructor.
 157
 158 //   if(fGeom)      delete fGeom;
 159 //   if(fCalibData) delete fCalibData;
 160 //   if(fCaloPed)   delete fCaloPed;
 161
 162 //   if (fDigitsArr) {
 163 //     fDigitsArr->Clear("C");
 164 //     delete fDigitsArr;
 165 //   }
 166 //   if (fRecPoints) {
 167 //     fRecPoints->Delete();
 168 //     delete fRecPoints;
 169 //  }
 170
 171   if(fClusterUnfolding) delete fClusterUnfolding;
 172 }
 173
 174 //____________________________________________________________________________
 175 Float_t  AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId)
 176 {
 177
 178   // Convert digitized amplitude into energy.
 179   // Calibration parameters are taken from calibration data base for raw data,
 180   // or from digitizer parameters for simulated data.
 181   if(fCalibData){
 182
 183     if (fGeom==0)
 184       AliFatal("Did not get geometry from EMCALLoader") ;
 185
 186     Int_t iSupMod = -1;
 187     Int_t nModule = -1;
 188     Int_t nIphi   = -1;
 189     Int_t nIeta   = -1;
 190     Int_t iphi    = -1;
 191     Int_t ieta    = -1;
 192
 193     Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
 194     if(!bCell) {
 195       fGeom->PrintGeometry();
 196       Error("Calibrate()"," Wrong cell id number : %i", absId);
 197       assert(0);
 198     }
 199
 200     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 201
 202     // Check if channel is bad (dead or hot), in this case return 0.
 203     // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
 204     // for the moment keep it here but remember to do the selection at the sdigitizer level
 205     // and remove it from here
 206     Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
 207     if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
 208                   AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
 209                   return 0;
 210     }
 211     //Check if time is too large or too small, indication of a noisy channel, remove in this case
 212     if(time > fTimeMax || time < fTimeMin) return 0;
 213
 214     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 215     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 216
 217     return -fADCpedestalECA + amp * fADCchannelECA ;
 218
 219   }
 220   else //Return energy with default parameters if calibration is not available
 221     return -fADCpedestalECA + amp * fADCchannelECA ;
 222
 223 }
 224
 225 //____________________________________________________________________________
 226 void AliEMCALClusterizer::GetCalibrationParameters()
 227 {
 228   // Set calibration parameters:
 229   // if calibration database exists, they are read from database,
 230   // otherwise, they are taken from digitizer.
 231   //
 232   // It is a user responsilibity to open CDB before reconstruction,
 233   // for example:
 234   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 235
 236   //Check if calibration is stored in data base
 237
 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 //____________________________________________________________________________
 251 void AliEMCALClusterizer::GetCaloCalibPedestal()
 252 {
 253   // Set calibration parameters:
 254   // if calibration database exists, they are read from database,
 255   // otherwise, they are taken from digitizer.
 256   //
 257   // It is a user responsilibity to open CDB before reconstruction,
 258   // for example:
 259   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 260
 261   //Check if calibration is stored in data base
 262
 263   if(!fCaloPed)
 264     {
 265       AliCDBEntry *entry = (AliCDBEntry*)
 266         AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 267       if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
 268     }
 269
 270   if(!fCaloPed)
 271     AliFatal("Pedestal info not found in CDB!");
 272
 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 //____________________________________________________________________________
 307 void AliEMCALClusterizer::InitParameters()
 308 {
 309   // Initializes the parameters for the Clusterizer
 310   fNumberOfECAClusters = 0 ;
 311   fCalibData           = 0 ;
 312   fCaloPed             = 0 ;
 313
 314   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 315   if(!recParam) {
 316     AliFatal("Reconstruction parameters for EMCAL not set!");
 317   } else {
 318     fECAClusteringThreshold = recParam->GetClusteringThreshold();
 319     fECAW0                  = recParam->GetW0();
 320     fMinECut                = recParam->GetMinECut();
 321     fToUnfold               = recParam->GetUnfold();
 322     if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");
 323     fECALocMaxCut           = recParam->GetLocMaxCut();
 324     fTimeCut                = recParam->GetTimeCut();
 325     fTimeMin                = recParam->GetTimeMin();
 326     fTimeMax                = recParam->GetTimeMax();
 327
 328     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",
 329                     fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
 330   }
 331
 332   if(fToUnfold){
 333     Int_t i=0;
 334     for (i = 0; i < 8; i++) {
 335       fSSPars[i] = recParam->GetSSPars(i);
 336     }//end of loop over parameters
 337     for (i = 0; i < 3; i++) {
 338       fPar5[i] = recParam->GetPar5(i);
 339       fPar6[i] = recParam->GetPar6(i);
 340     }//end of loop over parameters
 341
 342     fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6);
 343
 344     for (i = 0; i < 8; i++) {
 345       AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
 346     }
 347     for (i = 0; i < 3; i++) {
 348       AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
 349       AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
 350     }
 351
 352   }
 353
 354 }
 355
 356
 357 //____________________________________________________________________________
 358 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
 359 {
 360   // Print clusterizer parameters
 361
 362   TString message("\n") ;
 363
 364   if( strcmp(GetName(), "") !=0 ){
 365
 366     // Print parameters
 367
 368     TString taskName(Version()) ;
 369
 370     printf("--------------- ");
 371     printf("%s",taskName.Data()) ;
 372     printf(" ");
 373     printf("Clusterizing digits: ");
 374     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 375     printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 376     if(fToUnfold){
 377       printf("\nUnfolding on\n");
 378       printf("Unfolding parameters: fSSpars: \n");
 379       Int_t i=0;
 380       for (i = 0; i < 8; i++) {
 381         printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
 382       }
 383       printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
 384       for (i = 0; i < 3; i++) {
 385         printf("fPar5[%d] = %f \n", i, fPar5[i]);
 386         printf("fPar6[%d] = %f \n", i, fPar6[i]);
 387       }
 388     }
 389     else
 390       printf("\nUnfolding off\n");
 391
 392     printf("------------------------------------------------------------------");
 393   }
 394   else
 395     printf("AliEMCALClusterizer not initialized ") ;
 396 }
 397
 398 //____________________________________________________________________________
 399 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
 400 {
 401   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 402   if(strstr(option,"deb")) {
 403     printf("PrintRecPoints: Clusterization result:") ;
 404
 405     printf("           Found %d ECA Rec Points\n ",
 406            fRecPoints->GetEntriesFast()) ;
 407   }
 408
 409   if(strstr(option,"all")) {
 410     if(strstr(option,"deb")) {
 411       printf("\n-----------------------------------------------------------------------\n") ;
 412       printf("Clusters in ECAL section\n") ;
 413       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 414     }
 415     Int_t index;
 416     for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
 417       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
 418       if(rp){
 419         TVector3  globalpos;
 420         //rp->GetGlobalPosition(globalpos);
 421         TVector3  localpos;
 422         rp->GetLocalPosition(localpos);
 423         Float_t lambda[2];
 424         rp->GetElipsAxis(lambda);
 425
 426         Int_t nprimaries=0;
 427         Int_t * primaries = rp->GetPrimaries(nprimaries);
 428
 429         if(strstr(option,"deb"))
 430           printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 431                  rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 432                  globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 433                  rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 434         if(strstr(option,"deb")){
 435           for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 436             printf("%d ", primaries[iprimary] ) ;
 437           }
 438         }
 439       }
 440     }
 441
 442     if(strstr(option,"deb"))
 443       printf("\n-----------------------------------------------------------------------\n");
 444   }
 445 }
 446
 447 //___________________________________________________________________
 448 void  AliEMCALClusterizer::PrintRecoInfo()
 449 {
 450   //print reco version
 451   printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
 452
 453 }
 454
 455 //____________________________________________________________________________
 456 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
 457 {
 458   // Read the digits from the input tree
 459   TBranch *branch = digitsTree->GetBranch("EMCAL");
 460   if (!branch) {
 461     AliError("can't get the branch with the EMCAL digits !");
 462     return;
 463   }
 464   if (!fDigitsArr)
 465     fDigitsArr = new TClonesArray("AliEMCALDigit",100);
 466   branch->SetAddress(&fDigitsArr);
 467   branch->GetEntry(0);
 468 }
 469
 470 //____________________________________________________________________________
 471 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
 472 {
 473   // Read the digits from the input tree
 474   fTreeR = clustersTree;
 475
 476   AliDebug(9, "Making array for EMCAL clusters");
 477   fRecPoints = new TObjArray(100) ;
 478   Int_t split = 0;
 479   Int_t bufsize = 32000;
 480   fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
 481 }
 482
 483
 484