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