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