PWG4/PartCorrBase/AliCalorimeterUtils.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 /* $Id:  $ */
  16
  17 //_________________________________________________________________________
  18 // Class utility for Calorimeter specific selection methods                ///
  19 //
  20 //
  21 //
  22 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25
  26 // --- ROOT system ---
  27 #include "TGeoManager.h"
  28
  29 //---- ANALYSIS system ----
  30 #include "AliCalorimeterUtils.h"
  31 #include "AliESDEvent.h"
  32 #include "AliMCEvent.h"
  33 #include "AliStack.h"
  34 #include "AliAODPWG4Particle.h"
  35 #include "AliVCluster.h"
  36 #include "AliVCaloCells.h"
  37 #include "AliMixedEvent.h"
  38
  39 ClassImp(AliCalorimeterUtils)
  40
  41
  42 //____________________________________________
  43   AliCalorimeterUtils::AliCalorimeterUtils() :
  44     TObject(), fDebug(0),
  45     fEMCALGeoName("EMCAL_COMPLETEV1"),fPHOSGeoName("PHOSgeo"),
  46     fEMCALGeo(0x0),                   fPHOSGeo(0x0),
  47     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE),
  48     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
  49     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0),
  50     fNCellsFromPHOSBorder(0),
  51     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
  52     fRecalibration(kFALSE),           fPHOSRecalibrationFactors(),
  53     fEMCALRecoUtils(new AliEMCALRecoUtils),
  54     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
  55     fRecalculateMatching(kFALSE),
  56     fCutR(20),                        fCutZ(20),
  57     fCutEta(20),                      fCutPhi(20)
  58 {
  59   //Ctor
  60
  61   //Initialize parameters
  62   InitParameters();
  63   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
  64   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
  65
  66 }
  67
  68 //_________________________________________
  69 AliCalorimeterUtils::~AliCalorimeterUtils()
  70 {
  71   //Dtor
  72
  73   //if(fPHOSGeo)  delete fPHOSGeo  ;
  74   if(fEMCALGeo) delete fEMCALGeo ;
  75
  76   if(fPHOSBadChannelMap) {
  77     fPHOSBadChannelMap->Clear();
  78     delete  fPHOSBadChannelMap;
  79   }
  80
  81   if(fPHOSRecalibrationFactors) {
  82     fPHOSRecalibrationFactors->Clear();
  83     delete  fPHOSRecalibrationFactors;
  84   }
  85
  86   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
  87   if(fNMaskCellColumns) delete [] fMaskCellColumns;
  88
  89 }
  90
  91 //_____________________________________________________________________________________
  92 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
  93                                                     AliVCaloCells* cells,
  94                                                     AliVEvent * event, Int_t iev) const
  95 {
  96
  97         // Given the list of AbsId of the cluster, get the maximum cell and
  98         // check if there are fNCellsFromBorder from the calorimeter border
  99
 100   //If the distance to the border is 0 or negative just exit accept all clusters
 101         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
 102         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 103
 104   Int_t absIdMax        = -1;
 105         Float_t ampMax  = -1;
 106
 107   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 108   Int_t nMixedEvents = 0 ;
 109   Int_t * cellsCumul = NULL ;
 110   Int_t numberOfCells = 0 ;
 111   if (mixEvent){
 112     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 113     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 114       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 115       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 116     }
 117
 118     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 119       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 120       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 121     }
 122
 123     if(cellsCumul){
 124
 125       Int_t startCell = cellsCumul[iev] ;
 126       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 127       //Find cells with maximum amplitude
 128       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 129         Int_t absId = cluster->GetCellAbsId(i) ;
 130         for (Int_t j = startCell; j < endCell ;  j++) {
 131           Short_t cellNumber;
 132           Double_t amp ;
 133           Double_t time;
 134           cells->GetCell(j, cellNumber, amp, time) ;
 135           if (absId == cellNumber) {
 136             if(amp > ampMax){
 137               ampMax   = amp;
 138               absIdMax = absId;
 139             }
 140           }
 141         }
 142       }//loop on cluster cells
 143     }// cells cumul available
 144     else {
 145       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 146       abort();
 147     }
 148   } else {//Normal SE Events
 149     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 150       Int_t absId = cluster->GetCellAbsId(i) ;
 151       Float_t amp       = cells->GetCellAmplitude(absId);
 152       if(amp > ampMax){
 153         ampMax   = amp;
 154         absIdMax = absId;
 155       }
 156     }
 157   }
 158
 159         if(fDebug > 1)
 160                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 161            absIdMax, ampMax, cluster->E());
 162
 163         if(absIdMax==-1) return kFALSE;
 164
 165         //Check if the cell is close to the borders:
 166         Bool_t okrow = kFALSE;
 167         Bool_t okcol = kFALSE;
 168
 169         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 170
 171                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 172                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 173                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 174                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 175       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 176     }
 177
 178                 //Check rows/phi
 179     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 180                 if(iSM < 10){
 181                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 182     }
 183                 else{
 184                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
 185                 }
 186
 187                 //Check columns/eta
 188                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
 189                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 190                 }
 191                 else{
 192                         if(iSM%2==0){
 193                                 if(ieta >= nborder)     okcol = kTRUE;
 194                         }
 195                         else {
 196                                 if(ieta <  48-nborder)  okcol = kTRUE;
 197                         }
 198                 }//eta 0 not checked
 199                 if(fDebug > 1)
 200                 {
 201                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 202              nborder, ieta, iphi, iSM);
 203                         if (okcol && okrow ) printf(" YES \n");
 204                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 205                 }
 206         }//EMCAL
 207         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 208                 Int_t relId[4];
 209                 Int_t irow = -1, icol = -1;
 210                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 211                 irow = relId[2];
 212                 icol = relId[3];
 213                 //imod = relId[0]-1;
 214                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 215                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 216                 if(fDebug > 1)
 217                 {
 218                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 219              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 220                         if (okcol && okrow ) printf(" YES \n");
 221                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 222                 }
 223         }//PHOS
 224
 225         if (okcol && okrow) return kTRUE;
 226         else                return kFALSE;
 227
 228 }
 229
 230 //_________________________________________________________________________________________________________
 231 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
 232 {
 233         // Check that in the cluster cells, there is no bad channel of those stored
 234         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 235
 236         if (!fRemoveBadChannels) return kFALSE;
 237         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 238         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 239         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 240
 241         Int_t icol = -1;
 242         Int_t irow = -1;
 243         Int_t imod = -1;
 244         for(Int_t iCell = 0; iCell<nCells; iCell++){
 245
 246                 //Get the column and row
 247                 if(calorimeter == "EMCAL"){
 248       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 249                 }
 250                 else if(calorimeter=="PHOS"){
 251                         Int_t    relId[4];
 252                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 253                         irow = relId[2];
 254                         icol = relId[3];
 255                         imod = relId[0]-1;
 256                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 257                         if(GetPHOSChannelStatus(imod, icol, irow)) return kTRUE;
 258                 }
 259                 else return kFALSE;
 260
 261         }// cell cluster loop
 262
 263         return kFALSE;
 264
 265 }
 266
 267 //_______________________________________________________________
 268 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
 269 {
 270   // Correct cluster energy non linearity
 271
 272   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 273
 274 }
 275
 276 //___________________________________________________________________________________________________________________
 277 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, AliVCluster* clu, Float_t & clusterFraction) const
 278 {
 279
 280   //For a given CaloCluster gets the absId of the cell
 281   //with maximum energy deposit.
 282
 283   if( !clu || !cells ){
 284     AliInfo("Cluster or cells pointer is null!");
 285     return -1;
 286   }
 287
 288   Double_t eMax        =-1.;
 289   Double_t eTot        = 0.;
 290   Double_t eCell       =-1.;
 291   Float_t  fraction    = 1.;
 292   Float_t  recalFactor = 1.;
 293   Int_t    cellAbsId   =-1 , absId =-1 ;
 294   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
 295
 296   TString           calo = "EMCAL";
 297   if(clu->IsPHOS()) calo = "PHOS";
 298
 299   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 300
 301     cellAbsId = clu->GetCellAbsId(iDig);
 302
 303     fraction  = clu->GetCellAmplitudeFraction(iDig);
 304     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 305
 306     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
 307
 308     if(IsRecalibrationOn()) {
 309       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 310       else              recalFactor = GetPHOSChannelRecalibrationFactor(iSupMod,ieta,iphi);
 311     }
 312
 313     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 314
 315     if(eCell > eMax)  {
 316       eMax  = eCell;
 317       absId = cellAbsId;
 318     }
 319
 320     eTot+=eCell;
 321
 322   }// cell loop
 323
 324   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
 325
 326   return absId;
 327
 328 }
 329
 330 //_____________________________________________________________________________________________________
 331 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 332 {
 333         //Get the EMCAL/PHOS module number that corresponds to this particle
 334
 335         Int_t absId = -1;
 336         if(particle->GetDetector()=="EMCAL"){
 337                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 338                 if(fDebug > 2)
 339                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 340              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 341                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 342         }//EMCAL
 343         else if(particle->GetDetector()=="PHOS"){
 344     // In case we use the MC reader, the input are TParticles,
 345     // in this case use the corresponing method in PHOS Geometry to get the particle.
 346     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 347     {
 348       Int_t mod =-1;
 349       Double_t z = 0., x=0.;
 350       TParticle* primary = 0x0;
 351       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 352       if(stack) {
 353         primary = stack->Particle(particle->GetCaloLabel(0));
 354       }
 355       else {
 356         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 357       }
 358
 359       if(primary){
 360         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 361       }
 362       else{
 363         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 364       }
 365       return mod;
 366     }
 367     // Input are ESDs or AODs, get the PHOS module number like this.
 368     else{
 369       //FIXME
 370       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 371       //return GetModuleNumber(cluster);
 372       //MEFIX
 373       return -1;
 374     }
 375         }//PHOS
 376
 377         return -1;
 378 }
 379
 380 //_____________________________________________________________________
 381 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 382 {
 383         //Get the EMCAL/PHOS module number that corresponds to this cluster
 384         TLorentzVector lv;
 385         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
 386   if(!cluster){
 387     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 388     return -1;
 389   }
 390         cluster->GetMomentum(lv,v);
 391         Float_t phi = lv.Phi();
 392         if(phi < 0) phi+=TMath::TwoPi();
 393         Int_t absId = -1;
 394         if(cluster->IsEMCAL()){
 395                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
 396                 if(fDebug > 2)
 397                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 398              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 399                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 400         }//EMCAL
 401         else if(cluster->IsPHOS()) {
 402                 Int_t    relId[4];
 403                 if ( cluster->GetNCells() > 0) {
 404                         absId = cluster->GetCellAbsId(0);
 405                         if(fDebug > 2)
 406                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
 407                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
 408                 }
 409                 else return -1;
 410
 411                 if ( absId >= 0) {
 412                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 413                         if(fDebug > 2)
 414                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
 415                         return relId[0]-1;
 416                 }
 417                 else return -1;
 418         }//PHOS
 419
 420         return -1;
 421 }
 422
 423 //___________________________________________________________________________________________________
 424 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
 425                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
 426 {
 427         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 428         Int_t imod = -1;
 429         if ( absId >= 0) {
 430                 if(calo=="EMCAL"){
 431                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 432                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 433                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 434       if(imod < 0 || irow < 0 || icol < 0 ) {
 435         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
 436       }
 437
 438                         //RCU0
 439                         if (0<=irow&&irow<8) iRCU=0; // first cable row
 440                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
 441                         //second cable row
 442                         //RCU1
 443                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
 444                         //second cable row
 445                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
 446
 447                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 448                         if (iRCU<0) {
 449                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
 450                         }
 451
 452                         return imod ;
 453                 }//EMCAL
 454                 else{//PHOS
 455                         Int_t    relId[4];
 456                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 457                         irow = relId[2];
 458                         icol = relId[3];
 459                         imod = relId[0]-1;
 460                         iRCU= (Int_t)(relId[2]-1)/16 ;
 461                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
 462                         if (iRCU >= 4) {
 463                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
 464                         }
 465                         return imod;
 466                 }//PHOS
 467         }
 468
 469         return -1;
 470 }
 471
 472 //________________________________________
 473 void AliCalorimeterUtils::InitParameters()
 474 {
 475   //Initialize the parameters of the analysis.
 476   fEMCALGeoName         = "EMCAL_COMPLETEV1";
 477   fPHOSGeoName          = "PHOSgeo";
 478   fEMCALGeoMatrixSet    = kFALSE;
 479   fPHOSGeoMatrixSet     = kFALSE;
 480   fRemoveBadChannels    = kFALSE;
 481   fNCellsFromPHOSBorder = 0;
 482
 483   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
 484   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
 485   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
 486   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
 487   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 488   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 489
 490 }
 491
 492
 493 //_____________________________________________________
 494 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
 495 {
 496   //Init PHOS bad channels map
 497   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
 498   //In order to avoid rewriting the same histograms
 499   Bool_t oldStatus = TH1::AddDirectoryStatus();
 500   TH1::AddDirectory(kFALSE);
 501
 502   fPHOSBadChannelMap = new TObjArray(5);
 503   for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOSBadChannelMap_Mod%d",i),Form("PHOSBadChannelMap_Mod%d",i), 56, 0, 56, 64, 0, 64));
 504
 505   fPHOSBadChannelMap->SetOwner(kTRUE);
 506   fPHOSBadChannelMap->Compress();
 507
 508   //In order to avoid rewriting the same histograms
 509   TH1::AddDirectory(oldStatus);
 510 }
 511
 512 //______________________________________________________
 513 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
 514 {
 515         //Init EMCAL recalibration factors
 516         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
 517         //In order to avoid rewriting the same histograms
 518         Bool_t oldStatus = TH1::AddDirectoryStatus();
 519         TH1::AddDirectory(kFALSE);
 520
 521         fPHOSRecalibrationFactors = new TObjArray(5);
 522         for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 56, 0, 56, 64, 0, 64));
 523         //Init the histograms with 1
 524         for (Int_t m = 0; m < 5; m++) {
 525                 for (Int_t i = 0; i < 56; i++) {
 526                         for (Int_t j = 0; j < 64; j++) {
 527                                 SetPHOSChannelRecalibrationFactor(m,i,j,1.);
 528                         }
 529                 }
 530         }
 531         fPHOSRecalibrationFactors->SetOwner(kTRUE);
 532         fPHOSRecalibrationFactors->Compress();
 533
 534         //In order to avoid rewriting the same histograms
 535         TH1::AddDirectory(oldStatus);
 536 }
 537
 538
 539 //___________________________________________
 540 void AliCalorimeterUtils::InitEMCALGeometry()
 541 {
 542         //Initialize EMCAL geometry if it did not exist previously
 543         if (!fEMCALGeo){
 544                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 545                 if(fDebug > 0){
 546                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
 547                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 548                         printf("\n");
 549                 }
 550         }
 551 }
 552
 553 //__________________________________________
 554 void AliCalorimeterUtils::InitPHOSGeometry()
 555 {
 556         //Initialize PHOS geometry if it did not exist previously
 557         if (!fPHOSGeo){
 558                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
 559                 if(fDebug > 0){
 560                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
 561                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 562                         printf("\n");
 563                 }
 564         }
 565 }
 566
 567 //_________________________________________________________
 568 void AliCalorimeterUtils::Print(const Option_t * opt) const
 569 {
 570
 571   //Print some relevant parameters set for the analysis
 572   if(! opt)
 573     return;
 574
 575   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 576   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
 577   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
 578          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
 579   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
 580   printf("Recalibrate Clusters? %d\n",fRecalibration);
 581   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
 582   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 583   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
 584
 585   printf("    \n") ;
 586 }
 587
 588 //__________________________________________________________________
 589 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
 590                                              const Int_t ieta) const
 591 {
 592   //Check if cell is in one of the regions where we have significant amount
 593   //of material in front. Only EMCAL
 594
 595   Int_t icol = ieta;
 596   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 597
 598   if (fNMaskCellColumns && fMaskCellColumns) {
 599     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 600       if(icol==fMaskCellColumns[imask]) return kTRUE;
 601     }
 602   }
 603
 604   return kFALSE;
 605
 606 }
 607
 608
 609 //__________________________________________________________________________
 610 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
 611                                                       AliVCaloCells * cells)
 612 {
 613         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 614
 615   //Initialize some used variables
 616         Float_t energy   = 0;
 617         Int_t absId      = -1;
 618         Int_t icol = -1, irow = -1, iRCU = -1, module=1;
 619         Float_t factor = 1, frac = 0;
 620
 621         if(cells) {
 622
 623     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 624     UShort_t * index    = cluster->GetCellsAbsId() ;
 625     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 626     Int_t ncells     = cluster->GetNCells();
 627     TString calo     = "EMCAL";
 628     if(cluster->IsPHOS()) calo = "PHOS";
 629
 630     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 631     for(Int_t icell = 0; icell < ncells; icell++){
 632       absId = index[icell];
 633       frac =  fraction[icell];
 634       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
 635       module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
 636       if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow);
 637       else                  factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
 638       if(fDebug>2)
 639         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n",
 640                calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
 641
 642       energy += cells->GetCellAmplitude(absId)*factor*frac;
 643     }
 644
 645     if(fDebug>1)
 646       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
 647
 648         }// cells available
 649   else{
 650     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
 651   }
 652
 653         return energy;
 654 }
 655
 656 //________________________________________________________________________________
 657 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
 658 {
 659   //Set the calorimeters transformation matrices
 660
 661   //Get the EMCAL transformation geometry matrices from ESD
 662   if(!fEMCALGeoMatrixSet && fEMCALGeo){
 663     if(fLoadEMCALMatrices){
 664       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
 665       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 666         if(fEMCALMatrix[mod]){
 667           if(fDebug > 1)
 668             fEMCALMatrix[mod]->Print();
 669           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 670         }
 671       }//SM loop
 672       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 673
 674     }//Load matrices
 675     else if (!gGeoManager) {
 676
 677       if(fDebug > 1)
 678         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
 679       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 680         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 681           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 682           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
 683             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 684           }
 685         }// loop over super modules
 686         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 687
 688       }//ESD as input
 689       else {
 690         if(fDebug > 1)
 691           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 692       }//AOD as input
 693     }//Get matrix from data
 694     else if(gGeoManager){
 695       fEMCALGeoMatrixSet = kTRUE;
 696     }
 697   }//EMCAL geo && no geoManager
 698
 699         //Get the PHOS transformation geometry matrices from ESD
 700   if(!fPHOSGeoMatrixSet && fPHOSGeo){
 701     if(fLoadPHOSMatrices){
 702       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
 703       for(Int_t mod=0; mod < 5; mod++){
 704         if(fPHOSMatrix[mod]){
 705           if(fDebug > 1)
 706             fPHOSMatrix[mod]->Print();
 707           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
 708         }
 709       }//SM loop
 710       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
 711     }//Load matrices
 712     else if (!gGeoManager) {
 713       if(fDebug > 1)
 714         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
 715                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 716                                 for(Int_t mod=0; mod < 5; mod++){
 717                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
 718                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
 719                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
 720                                         }
 721                                 }// loop over modules
 722         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
 723                         }//ESD as input
 724                         else {
 725                                 if(fDebug > 1)
 726                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 727       }//AOD as input
 728     }// get matrix from data
 729     else if(gGeoManager){
 730       fPHOSGeoMatrixSet = kTRUE;
 731     }
 732         }//PHOS geo     and  geoManager was not set
 733 }
 734
 735 //__________________________________________________________________________________________
 736 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
 737 {
 738
 739   //Recalculate EMCAL cluster position
 740
 741   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
 742
 743 }
 744
 745 //________________________________________________________________________________
 746 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
 747                                                           TObjArray* clusterArray)
 748 {
 749   //Recalculate track matching
 750
 751   if (fRecalculateMatching) {
 752     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ;
 753     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
 754     //if(esdevent){
 755     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
 756     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
 757     //}
 758   }
 759 }