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
  16 //_________________________________________________________________________
  17 // Class utility for Calorimeter specific selection methods                ///
  18 //
  19 //
  20 //
  21 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 // --- ROOT system ---
  26 #include "TGeoManager.h"
  27
  28 // --- ANALYSIS system ---
  29 #include "AliCalorimeterUtils.h"
  30 #include "AliESDEvent.h"
  31 #include "AliMCEvent.h"
  32 #include "AliStack.h"
  33 #include "AliAODPWG4Particle.h"
  34 #include "AliVCluster.h"
  35 #include "AliVCaloCells.h"
  36 #include "AliMixedEvent.h"
  37
  38 // --- Detector ---
  39 #include "AliEMCALGeometry.h"
  40 #include "AliPHOSGeoUtils.h"
  41
  42 ClassImp(AliCalorimeterUtils)
  43
  44
  45 //____________________________________________
  46   AliCalorimeterUtils::AliCalorimeterUtils() :
  47     TObject(), fDebug(0),
  48     fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
  49     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     fLocMaxCutE(0),                   fLocMaxCutEDiff(0)
  63 {
  64   //Ctor
  65
  66   //Initialize parameters
  67   InitParameters();
  68   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
  69   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
  70
  71 }
  72
  73 //_________________________________________
  74 AliCalorimeterUtils::~AliCalorimeterUtils()
  75 {
  76   //Dtor
  77
  78   //if(fPHOSGeo)  delete fPHOSGeo  ;
  79   if(fEMCALGeo) delete fEMCALGeo ;
  80
  81   if(fPHOSBadChannelMap) {
  82     fPHOSBadChannelMap->Clear();
  83     delete  fPHOSBadChannelMap;
  84   }
  85
  86   if(fPHOSRecalibrationFactors) {
  87     fPHOSRecalibrationFactors->Clear();
  88     delete  fPHOSRecalibrationFactors;
  89   }
  90
  91   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
  92   if(fNMaskCellColumns) delete [] fMaskCellColumns;
  93
  94 }
  95
  96 //______________________________________________________________________________________
  97 Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
  98                                           const Int_t absId1, const Int_t absId2 ) const
  99 {
 100   // Tells if (true) or not (false) two cells are neighbours
 101   // A neighbour is defined as being two cells which share a side or corner
 102
 103   Bool_t areNeighbours = kFALSE ;
 104
 105   Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
 106   Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
 107
 108   Int_t rowdiff =  0, coldiff =  0;
 109
 110   Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
 111   Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
 112
 113   if(calo=="EMCAL" && nSupMod1!=nSupMod2)
 114   {
 115     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
 116     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 117     if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
 118     else           icol2+=AliEMCALGeoParams::fgkEMCALCols;
 119         }
 120
 121   rowdiff = TMath::Abs( irow1 - irow2 ) ;
 122   coldiff = TMath::Abs( icol1 - icol2 ) ;
 123
 124   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 125     areNeighbours = kTRUE ;
 126
 127   return areNeighbours;
 128
 129 }
 130
 131
 132 //_____________________________________________________________________________________
 133 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
 134                                                     AliVCaloCells* cells,
 135                                                     AliVEvent * event, Int_t iev) const
 136 {
 137
 138         // Given the list of AbsId of the cluster, get the maximum cell and
 139         // check if there are fNCellsFromBorder from the calorimeter border
 140
 141   //If the distance to the border is 0 or negative just exit accept all clusters
 142         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
 143         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 144
 145   Int_t absIdMax        = -1;
 146         Float_t ampMax  = -1;
 147
 148   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 149   Int_t nMixedEvents = 0 ;
 150   Int_t * cellsCumul = NULL ;
 151   Int_t numberOfCells = 0 ;
 152   if (mixEvent){
 153     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 154     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 155       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 156       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 157     }
 158
 159     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 160       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 161       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 162     }
 163
 164     if(cellsCumul){
 165
 166       Int_t startCell = cellsCumul[iev] ;
 167       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 168       //Find cells with maximum amplitude
 169       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 170         Int_t absId = cluster->GetCellAbsId(i) ;
 171         for (Int_t j = startCell; j < endCell ;  j++) {
 172           Short_t cellNumber;
 173           Double_t amp ;
 174           Double_t time;
 175           cells->GetCell(j, cellNumber, amp, time) ;
 176           if (absId == cellNumber) {
 177             if(amp > ampMax){
 178               ampMax   = amp;
 179               absIdMax = absId;
 180             }
 181           }
 182         }
 183       }//loop on cluster cells
 184     }// cells cumul available
 185     else {
 186       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 187       abort();
 188     }
 189   } else {//Normal SE Events
 190     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 191       Int_t absId = cluster->GetCellAbsId(i) ;
 192       Float_t amp       = cells->GetCellAmplitude(absId);
 193       if(amp > ampMax){
 194         ampMax   = amp;
 195         absIdMax = absId;
 196       }
 197     }
 198   }
 199
 200         if(fDebug > 1)
 201                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 202            absIdMax, ampMax, cluster->E());
 203
 204         if(absIdMax==-1) return kFALSE;
 205
 206         //Check if the cell is close to the borders:
 207         Bool_t okrow = kFALSE;
 208         Bool_t okcol = kFALSE;
 209
 210         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 211
 212                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 213                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 214                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 215                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 216       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 217     }
 218
 219                 //Check rows/phi
 220     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 221                 if(iSM < 10){
 222                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 223     }
 224                 else{
 225                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
 226                 }
 227
 228                 //Check columns/eta
 229                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
 230                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 231                 }
 232                 else{
 233                         if(iSM%2==0){
 234                                 if(ieta >= nborder)     okcol = kTRUE;
 235                         }
 236                         else {
 237                                 if(ieta <  48-nborder)  okcol = kTRUE;
 238                         }
 239                 }//eta 0 not checked
 240                 if(fDebug > 1)
 241                 {
 242                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 243              nborder, ieta, iphi, iSM);
 244                         if (okcol && okrow ) printf(" YES \n");
 245                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 246                 }
 247         }//EMCAL
 248         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 249                 Int_t relId[4];
 250                 Int_t irow = -1, icol = -1;
 251                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 252                 irow = relId[2];
 253                 icol = relId[3];
 254                 //imod = relId[0]-1;
 255                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 256                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 257                 if(fDebug > 1)
 258                 {
 259                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 260              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 261                         if (okcol && okrow ) printf(" YES \n");
 262                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 263                 }
 264         }//PHOS
 265
 266         if (okcol && okrow) return kTRUE;
 267         else                return kFALSE;
 268
 269 }
 270
 271 //_________________________________________________________________________________________________________
 272 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
 273 {
 274         // Check that in the cluster cells, there is no bad channel of those stored
 275         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 276
 277         if (!fRemoveBadChannels) return kFALSE;
 278         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 279         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 280         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 281
 282         Int_t icol = -1;
 283         Int_t irow = -1;
 284         Int_t imod = -1;
 285         for(Int_t iCell = 0; iCell<nCells; iCell++){
 286
 287                 //Get the column and row
 288                 if(calorimeter == "EMCAL"){
 289       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 290                 }
 291                 else if(calorimeter=="PHOS"){
 292                         Int_t    relId[4];
 293                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 294                         irow = relId[2];
 295                         icol = relId[3];
 296                         imod = relId[0]-1;
 297                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 298       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
 299                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
 300                 }
 301                 else return kFALSE;
 302
 303         }// cell cluster loop
 304
 305         return kFALSE;
 306
 307 }
 308
 309 //_______________________________________________________________
 310 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
 311 {
 312   // Correct cluster energy non linearity
 313
 314   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 315
 316 }
 317
 318 //________________________________________________________________________________________
 319 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
 320                                              Float_t & clusterFraction) const
 321 {
 322
 323   //For a given CaloCluster gets the absId of the cell
 324   //with maximum energy deposit.
 325
 326   if( !clu || !cells ){
 327     AliInfo("Cluster or cells pointer is null!");
 328     return -1;
 329   }
 330
 331   Double_t eMax        =-1.;
 332   Double_t eTot        = 0.;
 333   Double_t eCell       =-1.;
 334   Float_t  fraction    = 1.;
 335   Float_t  recalFactor = 1.;
 336   Int_t    cellAbsId   =-1 , absId =-1 ;
 337   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
 338
 339   TString           calo = "EMCAL";
 340   if(clu->IsPHOS()) calo = "PHOS";
 341
 342   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 343
 344     cellAbsId = clu->GetCellAbsId(iDig);
 345
 346     fraction  = clu->GetCellAmplitudeFraction(iDig);
 347     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 348
 349     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
 350
 351     if(IsRecalibrationOn()) {
 352       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 353       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
 354     }
 355
 356     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 357
 358     if(eCell > eMax)  {
 359       eMax  = eCell;
 360       absId = cellAbsId;
 361     }
 362
 363     eTot+=eCell;
 364
 365   }// cell loop
 366
 367   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
 368
 369   return absId;
 370
 371 }
 372
 373 //__________________________________________________________________________
 374 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
 375                                                  const AliVEvent* event,
 376                                                  const Int_t index) const
 377 {
 378   // Get the matched track given its index, usually just the first match
 379   // Since it is different for ESDs and AODs here it is a wrap method to do it
 380
 381   AliVTrack *track = 0;
 382
 383   // EMCAL case only when matching is recalculated
 384   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
 385   {
 386     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
 387     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
 388
 389     if(trackIndex < 0 )
 390     {
 391       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
 392     }
 393     else
 394     {
 395       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
 396     }
 397
 398     return track ;
 399
 400   }
 401
 402   // Normal case, get info from ESD or AOD
 403   // ESDs
 404   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 405   {
 406     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 407
 408     if(iESDtrack < 0 )
 409     {
 410       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
 411       return 0x0;
 412     }
 413
 414     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 415
 416   }
 417   else // AODs
 418   {
 419     if(cluster->GetNTracksMatched() > 0 )
 420       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
 421   }
 422
 423   return track ;
 424
 425 }
 426
 427 //_____________________________________________________________________________________________________
 428 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 429 {
 430         //Get the EMCAL/PHOS module number that corresponds to this particle
 431
 432         Int_t absId = -1;
 433         if(particle->GetDetector()=="EMCAL"){
 434                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 435                 if(fDebug > 2)
 436                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 437              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 438                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 439         }//EMCAL
 440         else if(particle->GetDetector()=="PHOS"){
 441     // In case we use the MC reader, the input are TParticles,
 442     // in this case use the corresponing method in PHOS Geometry to get the particle.
 443     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 444     {
 445       Int_t mod =-1;
 446       Double_t z = 0., x=0.;
 447       TParticle* primary = 0x0;
 448       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 449       if(stack) {
 450         primary = stack->Particle(particle->GetCaloLabel(0));
 451       }
 452       else {
 453         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 454       }
 455
 456       if(primary){
 457         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 458       }
 459       else{
 460         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 461       }
 462       return mod;
 463     }
 464     // Input are ESDs or AODs, get the PHOS module number like this.
 465     else{
 466       //FIXME
 467       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 468       //return GetModuleNumber(cluster);
 469       //MEFIX
 470       return -1;
 471     }
 472         }//PHOS
 473
 474         return -1;
 475 }
 476
 477 //_____________________________________________________________________
 478 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 479 {
 480         //Get the EMCAL/PHOS module number that corresponds to this cluster
 481         TLorentzVector lv;
 482         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
 483   if(!cluster){
 484     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 485     return -1;
 486   }
 487         cluster->GetMomentum(lv,v);
 488         Float_t phi = lv.Phi();
 489         if(phi < 0) phi+=TMath::TwoPi();
 490         Int_t absId = -1;
 491         if(cluster->IsEMCAL()){
 492                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
 493                 if(fDebug > 2)
 494                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 495              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 496                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 497         }//EMCAL
 498         else if(cluster->IsPHOS()) {
 499                 Int_t    relId[4];
 500                 if ( cluster->GetNCells() > 0) {
 501                         absId = cluster->GetCellAbsId(0);
 502                         if(fDebug > 2)
 503                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
 504                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
 505                 }
 506                 else return -1;
 507
 508                 if ( absId >= 0) {
 509                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 510                         if(fDebug > 2)
 511                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
 512                         return relId[0]-1;
 513                 }
 514                 else return -1;
 515         }//PHOS
 516
 517         return -1;
 518 }
 519
 520 //___________________________________________________________________________________________________
 521 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
 522                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
 523 {
 524         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 525         Int_t imod = -1;
 526         if ( absId >= 0) {
 527                 if(calo=="EMCAL"){
 528                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 529                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 530                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 531       if(imod < 0 || irow < 0 || icol < 0 ) {
 532         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
 533       }
 534
 535                         //RCU0
 536                         if (0<=irow&&irow<8) iRCU=0; // first cable row
 537                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
 538                         //second cable row
 539                         //RCU1
 540                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
 541                         //second cable row
 542                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
 543
 544                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 545                         if (iRCU<0) {
 546                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
 547                         }
 548
 549                         return imod ;
 550                 }//EMCAL
 551                 else{//PHOS
 552                         Int_t    relId[4];
 553                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 554                         irow = relId[2];
 555                         icol = relId[3];
 556                         imod = relId[0]-1;
 557                         iRCU= (Int_t)(relId[2]-1)/16 ;
 558                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
 559                         if (iRCU >= 4) {
 560                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
 561                         }
 562                         return imod;
 563                 }//PHOS
 564         }
 565
 566         return -1;
 567 }
 568
 569 //___________________________________________________________________________________________
 570 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
 571 {
 572   // Find local maxima in cluster
 573
 574   const Int_t   nc = cluster->GetNCells();
 575   Int_t   *absIdList = new Int_t  [nc];
 576   Float_t *maxEList  = new Float_t[nc];
 577
 578   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
 579
 580   delete [] absIdList;
 581   delete [] maxEList;
 582
 583   return nMax;
 584
 585 }
 586
 587 //___________________________________________________________________________________________
 588 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
 589                                                   Int_t *absIdList,     Float_t *maxEList)
 590 {
 591   // Find local maxima in cluster
 592
 593   Int_t iDigitN = 0 ;
 594   Int_t iDigit  = 0 ;
 595   Int_t absId1 = -1 ;
 596   Int_t absId2 = -1 ;
 597   const Int_t nCells = cluster->GetNCells();
 598
 599   TString calorimeter = "EMCAL";
 600   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
 601
 602   //printf("cluster : ncells %d \n",nCells);
 603   for(iDigit = 0; iDigit < nCells ; iDigit++)
 604   {
 605     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ;
 606      //Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 607      //RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
 608      //Int_t icol = -1, irow = -1, iRCU = -1;
 609      //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
 610      //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
 611   }
 612
 613   for(iDigit = 0 ; iDigit < nCells; iDigit++)
 614   {
 615     if(absIdList[iDigit]>=0)
 616     {
 617       absId1 = cluster->GetCellsAbsId()[iDigit];
 618
 619       Float_t en1 = cells->GetCellAmplitude(absId1);
 620       RecalibrateCellAmplitude(en1,calorimeter,absId1);
 621
 622       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
 623
 624       for(iDigitN = 0; iDigitN < nCells; iDigitN++)
 625       {
 626         absId2 = cluster->GetCellsAbsId()[iDigitN] ;
 627
 628         if(absId2==-1 || absId2==absId1) continue;
 629
 630         //printf("\t %d : absIDj %d\n",iDigitN, absId2);
 631
 632         Float_t en2 = cells->GetCellAmplitude(absId2);
 633         RecalibrateCellAmplitude(en2,calorimeter,absId2);
 634
 635         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
 636
 637         if ( AreNeighbours(calorimeter, absId1, absId2) )
 638         {
 639           // printf("\t \t Neighbours \n");
 640           if (en1 > en2 )
 641           {
 642             absIdList[iDigitN] = -1 ;
 643             //printf("\t \t indexN %d not local max\n",iDigitN);
 644             // but may be digit too is not local max ?
 645             if(en1 < en2 + fLocMaxCutEDiff) {
 646               //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
 647               absIdList[iDigit] = -1 ;
 648             }
 649           }
 650           else
 651           {
 652             absIdList[iDigit] = -1 ;
 653             //printf("\t \t index %d not local max\n",iDigitN);
 654             // but may be digitN too is not local max ?
 655             if(en1 > en2 - fLocMaxCutEDiff)
 656             {
 657               absIdList[iDigitN] = -1 ;
 658               //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
 659             }
 660           }
 661         } // if Are neighbours
 662         //else printf("\t \t NOT Neighbours \n");
 663       } // while digitN
 664     } // slot not empty
 665   } // while digit
 666
 667   iDigitN = 0 ;
 668   for(iDigit = 0; iDigit < nCells; iDigit++)
 669   {
 670     if(absIdList[iDigit]>=0 )
 671     {
 672       absIdList[iDigitN] = absIdList[iDigit] ;
 673       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 674       RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
 675       if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
 676       maxEList[iDigitN] = en ;
 677       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
 678       iDigitN++ ;
 679     }
 680   }
 681
 682   //printf("**********N maxima %d \n",iDigitN);
 683   //for(Int_t imax = 0; imax < iDigitN; imax++) printf("imax %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
 684
 685   return iDigitN ;
 686
 687 }
 688
 689 //________________________________________
 690 void AliCalorimeterUtils::InitParameters()
 691 {
 692   //Initialize the parameters of the analysis.
 693
 694   fEMCALGeoName         = "EMCAL_COMPLETEV1";
 695   fPHOSGeoName          = "PHOSgeo";
 696
 697   fEMCALGeoMatrixSet    = kFALSE;
 698   fPHOSGeoMatrixSet     = kFALSE;
 699
 700   fRemoveBadChannels    = kFALSE;
 701
 702   fNCellsFromPHOSBorder = 0;
 703
 704   fLocMaxCutE           = 0.1 ;
 705   fLocMaxCutEDiff       = 0.0 ;
 706
 707   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
 708   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
 709   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
 710   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
 711   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 712   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 713
 714 }
 715
 716
 717 //_____________________________________________________
 718 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
 719 {
 720   //Init PHOS bad channels map
 721   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
 722   //In order to avoid rewriting the same histograms
 723   Bool_t oldStatus = TH1::AddDirectoryStatus();
 724   TH1::AddDirectory(kFALSE);
 725
 726   fPHOSBadChannelMap = new TObjArray(5);
 727   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));
 728
 729   fPHOSBadChannelMap->SetOwner(kTRUE);
 730   fPHOSBadChannelMap->Compress();
 731
 732   //In order to avoid rewriting the same histograms
 733   TH1::AddDirectory(oldStatus);
 734 }
 735
 736 //______________________________________________________
 737 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
 738 {
 739         //Init EMCAL recalibration factors
 740         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
 741         //In order to avoid rewriting the same histograms
 742         Bool_t oldStatus = TH1::AddDirectoryStatus();
 743         TH1::AddDirectory(kFALSE);
 744
 745         fPHOSRecalibrationFactors = new TObjArray(5);
 746         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));
 747         //Init the histograms with 1
 748         for (Int_t m = 0; m < 5; m++) {
 749                 for (Int_t i = 0; i < 56; i++) {
 750                         for (Int_t j = 0; j < 64; j++) {
 751                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
 752                         }
 753                 }
 754         }
 755         fPHOSRecalibrationFactors->SetOwner(kTRUE);
 756         fPHOSRecalibrationFactors->Compress();
 757
 758         //In order to avoid rewriting the same histograms
 759         TH1::AddDirectory(oldStatus);
 760 }
 761
 762
 763 //___________________________________________
 764 void AliCalorimeterUtils::InitEMCALGeometry()
 765 {
 766         //Initialize EMCAL geometry if it did not exist previously
 767         if (!fEMCALGeo){
 768                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 769                 if(fDebug > 0){
 770                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
 771                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 772                         printf("\n");
 773                 }
 774         }
 775 }
 776
 777 //__________________________________________
 778 void AliCalorimeterUtils::InitPHOSGeometry()
 779 {
 780         //Initialize PHOS geometry if it did not exist previously
 781         if (!fPHOSGeo){
 782                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
 783                 if(fDebug > 0){
 784                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
 785                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 786                         printf("\n");
 787                 }
 788         }
 789 }
 790
 791 //_________________________________________________________
 792 void AliCalorimeterUtils::Print(const Option_t * opt) const
 793 {
 794
 795   //Print some relevant parameters set for the analysis
 796   if(! opt)
 797     return;
 798
 799   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 800   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
 801   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
 802          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
 803   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
 804   printf("Recalibrate Clusters? %d\n",fRecalibration);
 805   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
 806   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 807   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
 808
 809   printf("Loc. Max. E > %2.2f\n",       fLocMaxCutE);
 810   printf("Loc. Max. E Diff > %2.2f\n",  fLocMaxCutEDiff);
 811
 812   printf("    \n") ;
 813 }
 814
 815 //__________________________________________________________________
 816 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
 817                                              const Int_t ieta) const
 818 {
 819   //Check if cell is in one of the regions where we have significant amount
 820   //of material in front. Only EMCAL
 821
 822   Int_t icol = ieta;
 823   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 824
 825   if (fNMaskCellColumns && fMaskCellColumns) {
 826     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 827       if(icol==fMaskCellColumns[imask]) return kTRUE;
 828     }
 829   }
 830
 831   return kFALSE;
 832
 833 }
 834
 835 //__________________________________________________________________________________________
 836 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
 837                                                    const TString calo, const Int_t id) const
 838 {
 839   //Recaculate cell energy if recalibration factor
 840
 841   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
 842   Int_t nModule  = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
 843
 844   if (IsRecalibrationOn())
 845   {
 846     if(calo == "PHOS")
 847     {
 848       amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
 849     }
 850     else
 851     {
 852       amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
 853     }
 854   }
 855 }
 856
 857 //_________________________________________________________________________________
 858 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
 859                                               const TString calo,
 860                                               const Int_t id, const Int_t bc) const
 861 {
 862   // Recalculate time if time recalibration available for EMCAL
 863   // not ready for PHOS
 864
 865   if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
 866   {
 867     GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
 868   }
 869
 870 }
 871
 872
 873 //__________________________________________________________________________
 874 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
 875                                                       AliVCaloCells * cells)
 876 {
 877         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 878
 879   //Initialize some used variables
 880         Float_t frac  = 0., energy = 0.;
 881
 882         if(cells)
 883   {
 884     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 885
 886     UShort_t * index    = cluster->GetCellsAbsId() ;
 887     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 888
 889     Int_t ncells     = cluster->GetNCells();
 890
 891     TString calo     = "EMCAL";
 892     if(cluster->IsPHOS())
 893       calo = "PHOS";
 894
 895     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 896     for(Int_t icell = 0; icell < ncells; icell++){
 897
 898       Int_t absId = index[icell];
 899
 900       frac =  fraction[icell];
 901       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
 902
 903       Float_t amp = cells->GetCellAmplitude(absId);
 904       RecalibrateCellAmplitude(amp,calo, absId);
 905
 906       if(fDebug>2)
 907         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
 908                calo.Data(),frac,cells->GetCellAmplitude(absId));
 909
 910       energy += amp*frac;
 911     }
 912
 913     if(fDebug>1)
 914       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
 915
 916         }// cells available
 917   else
 918   {
 919     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
 920   }
 921
 922         return energy;
 923 }
 924
 925 //________________________________________________________________________________
 926 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
 927 {
 928   //Set the calorimeters transformation matrices
 929
 930   //Get the EMCAL transformation geometry matrices from ESD
 931   if(!fEMCALGeoMatrixSet && fEMCALGeo){
 932     if(fLoadEMCALMatrices){
 933       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
 934       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 935         if(fEMCALMatrix[mod]){
 936           if(fDebug > 1)
 937             fEMCALMatrix[mod]->Print();
 938           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 939         }
 940       }//SM loop
 941       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 942
 943     }//Load matrices
 944     else if (!gGeoManager) {
 945
 946       if(fDebug > 1)
 947         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
 948       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 949         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 950           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 951           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
 952             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 953           }
 954         }// loop over super modules
 955
 956         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 957
 958       }//ESD as input
 959       else {
 960         if(fDebug > 1)
 961           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 962       }//AOD as input
 963     }//Get matrix from data
 964     else if(gGeoManager){
 965       fEMCALGeoMatrixSet = kTRUE;
 966     }
 967   }//EMCAL geo && no geoManager
 968
 969         //Get the PHOS transformation geometry matrices from ESD
 970   if(!fPHOSGeoMatrixSet && fPHOSGeo){
 971
 972     if(fLoadPHOSMatrices){
 973       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
 974       for(Int_t mod = 0 ; mod < 5 ; mod++){
 975         if(fPHOSMatrix[mod]){
 976           if(fDebug > 1 )
 977             fPHOSMatrix[mod]->Print();
 978           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
 979         }
 980       }//SM loop
 981       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
 982
 983     }//Load matrices
 984     else if (!gGeoManager) {
 985       if(fDebug > 1)
 986         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
 987                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 988                                 for(Int_t mod = 0; mod < 5; mod++){
 989                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
 990                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
 991                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
 992                                         }
 993                                 }// loop over modules
 994         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
 995                         }//ESD as input
 996                         else {
 997                                 if(fDebug > 1)
 998                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 999       }//AOD as input
1000     }// get matrix from data
1001     else if(gGeoManager){
1002       fPHOSGeoMatrixSet = kTRUE;
1003     }
1004         }//PHOS geo     and  geoManager was not set
1005
1006 }
1007
1008 //__________________________________________________________________________________________
1009 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1010 {
1011
1012   //Recalculate EMCAL cluster position
1013
1014   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1015
1016 }
1017
1018 //________________________________________________________________________________
1019 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1020                                                           TObjArray* clusterArray)
1021 {
1022   //Recalculate track matching
1023
1024   if (fRecalculateMatching) {
1025     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ;
1026     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1027     //if(esdevent){
1028     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
1029     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1030     //}
1031   }
1032 }